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 global_4deg(pyOM):
  10    """ global 4 deg model with 15 levels
  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)    = (90,40,15)
  18      M.dt_mom    = 1800.0
  19      M.dt_tracer = 86400.0
  20 
  21      M.coord_degree     = 1
  22      M.enable_cyclic_x  = 1
  23 
  24      M.congr_epsilon = 1e-8
  25      M.congr_max_iterations = 20000
  26      M.enable_streamfunction = 1
  27 
  28      I=self.fortran.isoneutral_module   
  29      I.enable_neutral_diffusion = 1
  30      I.k_iso_0 = 1000.0
  31      I.k_iso_steep = 1000.0
  32      I.iso_dslope=4./1000.0
  33      I.iso_slopec=1./1000.0
  34      I.enable_skew_diffusion = 1
  35 
  36      M.enable_hor_friction  = 1; M.a_h = (4*M.degtom)**3*2e-11
  37      M.enable_hor_friction_cos_scaling = 1; M.hor_friction_cosPower=1
  38  
  39      M.enable_implicit_vert_friction = 1
  40      T=self.fortran.tke_module   
  41      T.enable_tke = 1
  42      T.c_k = 0.1
  43      T.c_eps = 0.7
  44      T.alpha_tke = 30.0
  45      T.mxl_min = 1e-8
  46      T.tke_mxl_choice = 2
  47      T.enable_tke_superbee_advection = 1
  48 
  49      E=self.fortran.eke_module   
  50      E.enable_eke = 1
  51      E.eke_k_max  = 1e4
  52      E.eke_c_k    = 0.4
  53      E.eke_c_eps  = 0.5
  54      E.eke_cross  = 2.
  55      E.eke_crhin  = 1.0
  56      E.eke_lmin   = 100.0
  57      E.enable_eke_superbee_advection = 1
  58 
  59      I=self.fortran.idemix_module   
  60      I.enable_idemix = 1
  61      I.enable_idemix_hor_diffusion = 1
  62      I.enable_eke_diss_surfbot = 1
  63      I.eke_diss_surfbot_frac = 0.2 # fraction which goes into bottom
  64      I.enable_idemix_superbee_advection = 1
  65 
  66      M.eq_of_state_type = 5
  67      return
  68 
  69 
  70    def set_grid(self):
  71        M=self.fortran.main_module   
  72        ddz = array([50.,70.,100.,140.,190.,240.,290.,340.,390.,440.,490.,540.,590.,640.,690.])
  73        M.dzt[:]=ddz[::-1]
  74        M.dxt[:] = 4.0
  75        M.dyt[:] = 4.0
  76        M.y_origin = -76.0
  77        M.x_origin = 4.0
  78        return
  79 
  80 
  81    def set_coriolis(self):
  82        M=self.fortran.main_module   
  83        for j in range( M.yt.shape[0] ): M.coriolis_t[:,j] = 2*M.omega*sin(M.yt[j]/180.*pi)
  84        return
  85 
  86    def set_topography(self):
  87          """ setup topography
  88          """
  89          M=self.fortran.main_module   
  90          M.kbot[:]=0
  91          bathy = reshape( fromfile('bathymetry.bin', dtype='>i'), (M.nx,M.ny), order='F' )
  92          lev_s = reshape( fromfile('lev_s.bin', dtype='>f', count=M.nx*M.ny*M.nz), (M.nx,M.ny,M.nz), order='F' )[:,:,::-1]
  93          for j in range(M.js_pe,M.je_pe+1):
  94            for i in range(M.is_pe,M.ie_pe+1):
  95               (ii,jj) = ( self.if2py(i), self.jf2py(j) )
  96               for k in range(M.nz-1,-1,-1): 
  97                   if lev_s[i-1,j-1,k] > 0.0: M.kbot[ii,jj]=k+1
  98               if bathy[i-1,j-1] == 0.0: M.kbot[ii,jj] =0
  99          M.kbot[ M.kbot == M.nz] = 0
 100          return
 101    
 102    def set_initial_conditions(self):
 103        """ setup initial conditions
 104        """
 105        M=self.fortran.main_module   
 106 
 107        self.taux = zeros( (M.i_blk+2*M.onx,M.j_blk+2*M.onx,12), 'd', order = 'F')
 108        self.tauy = zeros( (M.i_blk+2*M.onx,M.j_blk+2*M.onx,12), 'd', order = 'F')
 109        self.qnec = zeros( ( M.i_blk+2*M.onx, M.j_blk+2*M.onx, 12 ), 'd' , order = 'F')
 110        self.qnet = zeros( ( M.i_blk+2*M.onx, M.j_blk+2*M.onx, 12 ), 'd' , order = 'F')
 111        self.sst_clim = zeros( ( M.i_blk+2*M.onx, M.j_blk+2*M.onx, 12 ), 'd' , order = 'F')
 112        self.sss_clim = zeros( ( M.i_blk+2*M.onx, M.j_blk+2*M.onx, 12 ), 'd' , order = 'F')
 113        
 114        # initial conditions for T and S
 115        lev = reshape( fromfile('lev_t.bin', dtype='>f', count=M.nx*M.ny*M.nz), (M.nx,M.ny,M.nz), order='F' )
 116        lev = lev[M.is_pe-1:M.ie_pe,M.js_pe-1:M.je_pe,::-1]
 117        M.temp[M.onx:-M.onx,M.onx:-M.onx,:,M.tau-1] = lev*M.maskt[M.onx:-M.onx,M.onx:-M.onx,:]
 118        M.temp[...,M.taum1-1] = M.temp[...,M.tau-1]
 119        
 120        lev = reshape( fromfile('lev_s.bin', dtype='>f', count=M.nx*M.ny*M.nz), (M.nx,M.ny,M.nz), order='F' )
 121        lev = lev[M.is_pe-1:M.ie_pe,M.js_pe-1:M.je_pe,::-1]
 122        M.salt[M.onx:-M.onx,M.onx:-M.onx,:,M.tau-1] = lev*M.maskt[M.onx:-M.onx,M.onx:-M.onx,:]
 123        M.salt[...,M.taum1-1] = M.salt[...,M.tau-1]
 124 
 125        # use Trenberth wind stress from MITgcm instead of ECMWF (also contained in ecmwf_4deg.cdf)
 126        tx = reshape( fromfile('trenberth_taux.bin', dtype='>f'), (M.nx,M.ny,12), order='F' )/M.rho_0
 127        ty = reshape( fromfile('trenberth_tauy.bin', dtype='>f'), (M.nx,M.ny,12), order='F' )/M.rho_0
 128        self.taux[M.onx:-M.onx,M.onx:-M.onx,:] = tx[M.is_pe-1:M.ie_pe,M.js_pe-1:M.je_pe,:]
 129        self.tauy[M.onx:-M.onx,M.onx:-M.onx,:] = ty[M.is_pe-1:M.ie_pe,M.js_pe-1:M.je_pe,:]
 130 
 131        # heat flux 
 132        fid = NF('ecmwf_4deg_monthly.cdf','r')
 133        for k in range(12):
 134            self.qnec[M.onx:-M.onx,M.onx:-M.onx,k] = fid.variables['Q3'][k,M.js_pe-1:M.je_pe,M.is_pe-1:M.ie_pe].transpose()
 135        self.qnec[ self.qnec  <= -1e10 ] = 0.0
 136 
 137        q = reshape( fromfile('ncep_qnet.bin', dtype='>f'), (M.nx,M.ny,12), order='F' )
 138        self.qnet[M.onx:-M.onx,M.onx:-M.onx,:] = - q[M.is_pe-1:M.ie_pe,M.js_pe-1:M.je_pe,:]
 139        self.qnet[ self.qnet  <= -1e10 ] = 0.0
 140 
 141        fxa = zeros( (1,), 'd', order = 'F')
 142        fxb = zeros( (1,), 'd', order = 'F')
 143        for n in range(12):
 144            fxa = fxa + sum( self.qnet[M.onx:-M.onx,M.onx:-M.onx,n]*M.area_t[M.onx:-M.onx,M.onx:-M.onx] )
 145            fxb = fxb + sum( M.area_t[M.onx:-M.onx,M.onx:-M.onx] )
 146        self.fortran.global_sum(fxa)
 147        self.fortran.global_sum(fxb)
 148        fxa = fxa[0]/fxb[0]
 149        if M.my_pe==0: print  ' removing an annual mean heat flux imbalance of %e W/m^2'% fxa
 150        for n in range(12):
 151          self.qnet[:,:,n] = (self.qnet[:,:,n] - fxa)*M.maskt[:,:,-1]
 152 
 153        # SST and SSS
 154        sst = reshape( fromfile('lev_sst.bin', dtype='>f'), (M.nx,M.ny,12), order='F' )
 155        sss = reshape( fromfile('lev_sss.bin', dtype='>f'), (M.nx,M.ny,12), order='F' )
 156        self.sst_clim[M.onx:-M.onx,M.onx:-M.onx,:] = sst[M.is_pe-1:M.ie_pe,M.js_pe-1:M.je_pe,:]
 157        self.sss_clim[M.onx:-M.onx,M.onx:-M.onx,:] = sss[M.is_pe-1:M.ie_pe,M.js_pe-1:M.je_pe,:]
 158 
 159        I=self.fortran.idemix_module   
 160        if I.enable_idemix:
 161          f = reshape( fromfile('tidal_energy.bin', dtype='>f'), (M.nx,M.ny), order='F' )/M.rho_0
 162          I.forc_iw_bottom[ M.onx:-M.onx,M.onx:-M.onx]  = f[M.is_pe-1:M.ie_pe,M.js_pe-1:M.je_pe]
 163          f = reshape( fromfile('wind_energy.bin', dtype='>f'), (M.nx,M.ny), order='F' )/M.rho_0*0.2
 164          I.forc_iw_surface[M.onx:-M.onx,M.onx:-M.onx] = f[M.is_pe-1:M.ie_pe,M.js_pe-1:M.je_pe]
 165        return
 166 
 167    def get_periodic_interval(self,currentTime,cycleLength,recSpacing,nbRec):
 168        """  interpolation routine taken from mitgcm
 169        """
 170        locTime = currentTime - recSpacing*0.5 + cycleLength*( 2 - round(currentTime/cycleLength) )
 171        tmpTime = locTime % cycleLength 
 172        tRec1 = 1 + int( tmpTime/recSpacing )
 173        tRec2 = 1 + tRec1 % int(nbRec) 
 174        wght2 = ( tmpTime - recSpacing*(tRec1 - 1) )/recSpacing
 175        wght1 = 1.0 - wght2
 176        return (wght1,wght2,tRec1,tRec2)
 177    
 178 
 179    def set_forcing(self):
 180      M=self.fortran.main_module   
 181      
 182      (f1,f2,n1,n2) = self.get_periodic_interval((M.itt-1)*M.dt_tracer,365*86400.0,365*86400./12.,12)
 183 
 184      # wind stress
 185      M.surface_taux[:]=(f1*self.taux[:,:,n1-1] + f2*self.taux[:,:,n2-1])
 186      M.surface_tauy[:]=(f1*self.tauy[:,:,n1-1] + f2*self.tauy[:,:,n2-1])
 187 
 188      # tke flux
 189      T=self.fortran.tke_module   
 190      if T.enable_tke:
 191        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  \
 192                                             +(0.5*(M.surface_tauy[1:-1,1:-1]+M.surface_tauy[1:-1,:-2]) )**2 )**(3./2.) 
 193      # heat flux : W/m^2 K kg/J m^3/kg = K m/s 
 194      cp_0 = 3991.86795711963
 195      sst  =  f1*self.sst_clim[:,:,n1-1]+f2*self.sst_clim[:,:,n2-1] 
 196      qnec =  f1*self.qnec[:,:,n1-1]    +f2*self.qnec[:,:,n2-1]
 197      qnet =  f1*self.qnet[:,:,n1-1]    +f2*self.qnet[:,:,n2-1]
 198      M.forc_temp_surface[:] =(qnet+ qnec*(sst-M.temp[:,:,-1,M.tau-1]) )*M.maskt[:,:,-1]/cp_0/M.rho_0
 199 
 200      # salinity restoring
 201      t_rest= 30*86400.0
 202      sss  =  f1*self.sss_clim[:,:,n1-1]+f2*self.sss_clim[:,:,n2-1] 
 203      M.forc_salt_surface[:] =M.dzt[-1]/t_rest*(sss-M.salt[:,:,-1,M.tau-1])*M.maskt[:,:,-1]
 204 
 205      # apply simple ice mask                     
 206      n=nonzero( logical_and( M.temp[:,:,-1,M.tau-1]*M.maskt[:,:,-1] <= -1.8 , M.forc_temp_surface <= 0.0 ) )                            
 207      M.forc_temp_surface[n]=0.0
 208      M.forc_salt_surface[n]=0.0
 209 
 210      if M.enable_tempsalt_sources:
 211         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] )
 212         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] )
 213      return
 214 
 215    def set_diagnostics(self):
 216      M=self.fortran.main_module   
 217      self.register_average(name='temp',long_name='Temperature',         units = 'deg C' , grid = 'TTT', var = M.temp)
 218      self.register_average(name='salt',long_name='Salinity',            units = 'g/kg' ,  grid = 'TTT', var = M.salt)
 219      self.register_average(name='u',   long_name='Zonal velocity',      units = 'm/s' ,   grid = 'UTT', var = M.u)
 220      self.register_average(name='v',   long_name='Meridional velocity', units = 'm/s' ,   grid = 'TUT', var = M.v)
 221      self.register_average(name='w',   long_name='Vertical velocity',   units = 'm/s' ,   grid = 'TTU', var = M.w)
 222      self.register_average(name='taux',long_name='wind stress',         units = 'm^2/s' , grid = 'UT',  var = M.surface_taux)
 223      self.register_average(name='tauy',long_name='wind stress',         units = 'm^2/s' , grid = 'TU',  var = M.surface_tauy)
 224      self.register_average(name='psi' ,long_name='Streamfunction',      units = 'm^3/s' , grid = 'UU',  var = M.psi)
 225      return
 226 
 227 
 228    def user_defined_signal(self):
 229        """ this routine must be called by all processors
 230        """
 231        M=self.fortran.main_module  
 232        a = zeros( (M.nx,M.ny), 'd', order = 'F')
 233        a[M.is_pe-1:M.ie_pe,0] = M.xt[2:-2]
 234        self.fortran.pe0_recv_2d(a)
 235        self.xt_gl = a[:,0].copy()
 236        
 237        a[0,M.js_pe-1:M.je_pe] = M.yt[2:-2]
 238        self.fortran.pe0_recv_2d(a)
 239        self.yt_gl = a[0,:].copy()
 240        
 241        self.psi_gl = zeros( (M.nx,M.ny), 'd', order = 'F')
 242        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) 
 243        self.fortran.pe0_recv_2d(self.psi_gl)
 244        
 245        self.temp_gl = zeros( (M.nx,M.ny,M.nz), 'd', order = 'F')
 246        for k in range(M.nz):
 247          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) 
 248          self.fortran.pe0_recv_2d(a)
 249          self.temp_gl[:,:,k]=a.copy()
 250 
 251        self.salt_gl = zeros( (M.nx,M.ny,M.nz), 'd', order = 'F')
 252        for k in range(M.nz):
 253          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) 
 254          self.fortran.pe0_recv_2d(a)
 255          self.salt_gl[:,:,k]=a.copy()
 256 
 257        #self.kappa_gl = zeros( (M.nx,M.ny,M.nz), 'd', order = 'F')
 258        #for k in range(M.nz):
 259        #  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) 
 260        #  self.fortran.pe0_recv_2d(a)
 261        #  self.kappa_gl[:,:,k]=a.copy()
 262 
 263        return
 264 
 265  
 266 
 267    def make_plot(self):
 268        M=self.fortran.main_module         # fortran module with model variables
 269        self.figure.clf()
 270        
 271        self.set_signal('user_defined') # following routine is called by all PEs
 272        self.user_defined_signal()
 273        
 274        ax=self.figure.add_subplot(221)
 275        co=ax.contourf(self.xt_gl,self.yt_gl,self.psi_gl.transpose()*1e-6)
 276        self.figure.colorbar(co)
 277        ax.set_title('Streamfunction [Sv]')
 278        #ax.set_xlabel('Longitude [deg E]')
 279        ax.set_ylabel('Latitude [deg N]')
 280        ax.axis('tight')
 281        
 282        ax=self.figure.add_subplot(222)
 283        co=ax.contourf(self.xt_gl,self.yt_gl,self.temp_gl[:,:,-1].transpose())
 284        self.figure.colorbar(co)
 285        ax.set_title('Sea Surface Temperature [deg C]')
 286        #ax.set_xlabel('Longitude [deg E]')
 287        ax.set_ylabel('Latitude [deg N]')
 288        ax.axis('tight')
 289        
 290        ax=self.figure.add_subplot(223)
 291        co=ax.contourf(self.yt_gl,M.zt,self.temp_gl[M.nx/2-1,:,:].transpose())
 292        self.figure.colorbar(co)
 293        ax.set_title('Temperature [deg C]')
 294        ax.set_ylabel('z [m]')
 295        ax.axis('tight')
 296        
 297        ax=self.figure.add_subplot(224)
 298        co=ax.contourf(self.yt_gl,M.zt,self.salt_gl[M.nx/2-1,:,:].transpose())
 299        self.figure.colorbar(co)
 300        ax.set_title('Salinity [g/kg]')
 301        ax.set_ylabel('z [m]')
 302        ax.axis('tight')
 303        
 304        #ax=self.figure.add_subplot(224)
 305        #try:
 306        # co=ax.contourf(self.yt_gl,M.zw,log10(self.kappa_gl[M.nx/2-1,:,:].transpose()) )
 307        #except:
 308        # pass
 309        #self.figure.colorbar(co)
 310        #ax.set_title('Diffusivity')
 311        #ax.set_xlabel('Latitude [deg N]')
 312        #ax.set_ylabel('z [m]')
 313        #ax.axis('tight')
 314 
 315        
 316        return
 317 
 318 if __name__ == "__main__": global_4deg().run(snapint= 86400.0*25, runlen = 86400.*365*10)