Advertisement
Guest User

Untitled

a guest
Jan 23rd, 2017
83
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 5.14 KB | None | 0 0
  1. '''
  2. Compute and plot moc, as a post-processing computation
  3. August 2016, Mark Petersen, LANL
  4. '''
  5. from netCDF4 import Dataset
  6. import matplotlib.pyplot as plt
  7. import numpy as np
  8.  
  9. ########################################################################
  10. #
  11. # Set directory, file, and variable names
  12. #
  13. ########################################################################
  14. fileName = 'mpaso.hist.am.timeSeriesStats.yr100-119.avg.nc'
  15. fileNameMesh= 'mpaso.rst.0002-01-01_00000.nc'
  16. fileNameMask='EC60to30v3_SingleRegionAtlanticWTransportTransects_masks.nc'
  17. contourLines = np.arange(-8,20.1,2)
  18.  
  19. timeIndex = 0
  20. #horVelName = 'timeMonthly_avg_normalVelocity'
  21. #vertVelName = 'timeMonthly_avg_vertVelocityTop'
  22. horVelName = 'time_avg_normalVelocity'
  23. vertVelName = 'time_avg_vertVelocityTop'
  24.  
  25. # This is specific to Atlantic MOC:
  26. mocLat = np.arange(-34.5,70.1,0.5)
  27. #contourLines = np.arange(-6,18.1,2)
  28.  
  29. ########################################################################
  30. #
  31. # Open files, load variables
  32. #
  33. ########################################################################
  34. print('** load variables: \n')
  35.  
  36. print 'Reading data from: ' + fileName
  37. ncfile = Dataset(fileName,'r')
  38. nCells = ncfile.dimensions['nCells'].size
  39. nVertLevels = ncfile.dimensions['nVertLevels'].size
  40.  
  41. horVel = ncfile.variables[horVelName][:,:,:]
  42. vertVel = ncfile.variables[vertVelName][:,:,:]
  43.  
  44. # delete:
  45. #print np.shape(vertVel)
  46. #print np.size(vertVel)
  47.  
  48. print 'Reading data from: ' + fileNameMesh
  49. ncfileMesh = Dataset(fileNameMesh,'r')
  50. dvEdge = ncfileMesh.variables['dvEdge'][:]
  51. areaCell = ncfileMesh.variables['areaCell'][:]
  52. latCell = ncfileMesh.variables['latCell'][:]*180./np.pi
  53. refBottomDepth = ncfileMesh.variables['refBottomDepth'][:]
  54.  
  55. refLayerThickness = np.zeros(nVertLevels)
  56. refLayerThickness[0] = refBottomDepth[0]
  57. refLayerThickness[1:nVertLevels] = refBottomDepth[1:nVertLevels]-refBottomDepth[0:nVertLevels-1]
  58.  
  59. refTopDepth = np.zeros(nVertLevels+1)
  60. refTopDepth[1:nVertLevels+1] = refBottomDepth[0:nVertLevels]
  61.  
  62. print 'Reading transects from: ' + fileNameMask
  63. ncfileMask = Dataset(fileNameMask,'r')
  64. nTransects = ncfileMask.dimensions['nTransects'].size
  65. maxEdgesInTransect = ncfileMask.dimensions['maxEdgesInTransect'].size
  66.  
  67. transectEdgeMaskSigns = ncfileMask.variables['transectEdgeMaskSigns'][:,:]
  68. transectEdgeGlobalIDs = ncfileMask.variables['transectEdgeGlobalIDs'][:,:]
  69.  
  70. print 'Reading regionss from: ' + fileNameMask
  71. ncfileMask = Dataset(fileNameMask,'r')
  72. nRegions = ncfileMask.dimensions['nRegions'].size
  73. regionCellMasks = ncfileMask.variables['regionCellMasks']
  74.  
  75. ########################################################################
  76. #
  77. # Compute transport through transects
  78. #
  79. ########################################################################
  80. print('** compute transport: \n')
  81.  
  82. m3ps_to_Sv = 1e-6; # m^3/sec flux to Sverdrups
  83.  
  84. # the volume transport
  85. transport = np.zeros(nTransects);
  86. transportZ = np.zeros([nVertLevels,nTransects]);
  87. transportZEdge = np.zeros([nVertLevels,maxEdgesInTransect,nTransects]);
  88.  
  89. for iTransect in range(nTransects):
  90. for i in range(maxEdgesInTransect):
  91. if transectEdgeGlobalIDs[iTransect, i]==0:
  92. break
  93. iEdge = transectEdgeGlobalIDs[iTransect, i] - 1 # subtract 1 because of python 0-indexing
  94. for k in range(nVertLevels):
  95. transportZEdge[k,i,iTransect] = ( transectEdgeMaskSigns[iEdge,iTransect]
  96. * horVel[timeIndex,iEdge,k] * dvEdge[iEdge] * refLayerThickness[k]*m3ps_to_Sv )
  97. transportZ[k,iTransect] = transportZ[k,iTransect] + transportZEdge[k,i,iTransect]
  98. transport[iTransect] = transport[iTransect] + transportZEdge[k,i,iTransect]
  99.  
  100. ########################################################################
  101. #
  102. # Compute MOC
  103. #
  104. ########################################################################
  105. print('** compute moc: \n')
  106.  
  107. nLat = np.size(mocLat)
  108. mocTop = np.zeros((nLat,nVertLevels+1))
  109.  
  110. for iRegion in range(nRegions):
  111. # assume transects and regions have the same index ordering:
  112. iTransect = iRegion
  113.  
  114. for k in range(1,nVertLevels+1):
  115. mocTop[0,k] = mocTop[0,k-1] + transportZ[k-1,iTransect]/m3ps_to_Sv
  116.  
  117. for iLat in range(1,nLat):
  118. ind =np.logical_and(np.logical_and(regionCellMasks[:,iRegion]==1,latCell >= mocLat[iLat-1] ), latCell < mocLat[iLat])
  119.  
  120. for k in range(1,nVertLevels+1):
  121. mocTop[iLat,k] = mocTop[iLat-1,k] + np.sum(vertVel[timeIndex,ind,k]*areaCell[ind])
  122.  
  123. # convert m^3/s to Sverdrup
  124. mocTop = mocTop * m3ps_to_Sv
  125.  
  126. ########################################################################
  127. #
  128. # Plot MOC
  129. #
  130. ########################################################################
  131. print('** plot moc: \n')
  132.  
  133. plt.clf()
  134. plt.figure(figsize=(8,4.5))
  135. X, Y = np.meshgrid(mocLat,refTopDepth) # change to zMid
  136. contourSet = plt.contourf(X,Y,mocTop.T,levels=contourLines)
  137. contourSet2 = plt.contour(X,Y,mocTop.T,levels=contourLines[::2],colors='k',hold='on')
  138. plt.xlabel('latitude')
  139. plt.ylabel('depth, m')
  140. plt.title('MPAS-Ocean Atlantic MOC, Sv')
  141. plt.gca().invert_yaxis()
  142. #plt.clabel(contourSet, fontsize=10)
  143. # Make a colorbar for the ContourSet returned by the contourf call.
  144. cbar = plt.colorbar(contourSet)
  145. cbar.add_lines(contourSet2)
  146. imageFileName = 'moc'
  147. plt.savefig(imageFileName+'.png')
  148.  
  149. ncfile.close()
  150. ncfileMask.close()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement