<jupyter_start><jupyter_text>Model sensitivity to different values of horizontal and vertical prescribed diffusivity coefficients using GMREDI. All horizontal diffusivites have to be set to 0. KzTr was changed.  Sensitivity to Explicit Diffusivity Coefficients============================This ipython notebook investigates the model sensitivity to different values of vertical prescribed diffusivity coefficients using the GMREDI package. Experiment 360x360x90_8Tr_LinProfiles_BarkleyHyd_GMREDI has upwelling body forcing and a linear vertical gradient in the initial tracer concentration. Run01 explores values of vertical diffusivity (1E-7 to 10) keeping horizontal diffusivity fixed (0.0). OOPS! It only has 3 tracers because I didn't change the number of tracers flag in data.ptracers.<jupyter_code>#KRM
import os

import numpy as np

import matplotlib.pyplot as plt

from scipy.stats import nanmean

from math import *

%matplotlib inline

import numpy as np

import matplotlib.colors as mcolors

from netCDF4 import Dataset

import pylab as pl

from decimal import getcontext, Decimal, Inexact

import statsmodels.api as sm

import pandas as pd
#lib_path = os.path.abspath('../PythonModulesMITgcm') # Add absolute path to my python scripts
#sys.path.append(lib_path)

import NumDiffTools as dif<jupyter_output><empty_output><jupyter_text>Important input parameters:------------------------<jupyter_code>nx = 360  # number of x points in domain
ny = 360  # number of y points in domain
nz = 90  # number of z levels in domain
nt = 21   # number of output files (n time) 
delt= 12*3600.0 # time interval between model output<jupyter_output><empty_output><jupyter_text>Load output filenames:<jupyter_code>### Load output filenames ###

filename = '/ocean/kramosmu/MITgcm/CanyonUpwelling/360x360x90_8Tr_LinProfiles_BarkleyHyd_GMREDI/run01/state.0000000000.glob.nc'
StateOut = Dataset(filename)
#for dimobj in StateOut.variables.values():
#    print dimobj

filename2 = '/ocean/kramosmu/MITgcm/CanyonUpwelling/360x360x90_8Tr_LinProfiles_BarkleyHyd_GMREDI/run01/grid.glob.nc'
GridOut = Dataset(filename2)
for dimobj in GridOut.variables.values():
    print dimobj

filename3 = '/ocean/kramosmu/MITgcm/CanyonUpwelling/360x360x90_8Tr_LinProfiles_BarkleyHyd_GMREDI/run01/ptracers.0000000000.glob.nc'
PtracersOut = Dataset(filename3)
#for dimobj in PtracersOut.variables.values():
#    print dimobj<jupyter_output><type 'netCDF4.Variable'>
float64 Depth(Y, X)
    description: fluid thickness in r coordinates (at rest)
unlimited dimensions: 
current shape = (360, 360)

<type 'netCDF4.Variable'>
float64 X(X)
    long_name: X-coordinate of cell center
    units: meters
unlimited dimensions: 
current shape = (360,)

<type 'netCDF4.Variable'>
float64 Y(Y)
    long_name: Y-Coordinate of cell center
    units: meters
unlimited dimensions: 
current shape = (360,)

<type 'netCDF4.Variable'>
float64 Xp1(Xp1)
    long_name: X-Coordinate of cell corner
    units: meters
unlimited dimensions: 
current shape = (360,)

<type 'netCDF4.Variable'>
float64 dxC(Y, Xp1)
    description: x cell center separation
unlimited dimensions: 
current shape = (360, 360)

<type 'netCDF4.Variable'>
float64 dxF(Y, X)
    description: x cell face separation
unlimited dimensions: 
current shape = (360, 360)

<type 'netCDF4.Variable'>
float64 Yp1(Yp1)
    long_name: Y-Coordinate of cell corner
    units: meters
unlimited dimensions[...]<jupyter_text>Vertical diffusivity Sensitivity--------------------Let's calculate the horizontal and vertical diffusivity for each of the 3 tracers in run01.<jupyter_code>tracer1 = dif.getField(filename3,'Tr1')

tracer1msk = np.ma.masked_values(tracer1, 0)

MASK = np.ma.getmask(tracer1msk)
   
    
tracer2 = dif.getField(filename3,'Tr2')

tracer2msk = np.ma.masked_values(tracer2, 0)


tracer3 = dif.getField(filename3,'Tr3')

tracer3msk = np.ma.masked_values(tracer3, 0)<jupyter_output><empty_output><jupyter_text>Calculate the volume of the domain<jupyter_code>DomVol = 0.0
    
hFacC = dif.getField(filename2,'HFacC')
drF   = dif.getField(filename2,'drF')
rAC   = dif.getField(filename2,'rA')
    
for jj in np.arange(ny):
    for kk in np.arange(nz):
      DomVol = DomVol + hFacC[kk,jj,:]*drF[kk]*rAC[jj,:]
    #print DomVol
  #print('\n The total volume is %e m^3 \n' %DomVolume)
   
DomVolume=sum(DomVol)
print DomVolume  
#vol = dif.CalcDomVolume(filename2, nx, ny, nz)

#DomVolume = vol[0]
#hFacC = vol[1]
#drF = vol[2]
#rAC = vol[3]

#print(DomVolume)<jupyter_output>5.96159528708e+12<jupyter_text>Tracer 1<jupyter_code>q2mean1 = dif.CalcVariance(nt, nz, ny, nx, DomVolume, hFacC, drF, rAC, tracer1)
print('done')
dqdt1 = dif.CalcTimeDer(q2mean1, nt, delt)
print('done')
qmeanDh1 = dif.CalcAvgHorGrad(filename2, nt,nz,ny,nx,tracer1,DomVolume,hFacC, drF, rAC)
print('done')
qmeanDz1=dif.CalcAvgVerGrad(filename2, nt,nz,ny,nx,tracer1,DomVolume,hFacC, drF, rAC)
print('done')

esth1=dif.FitDiffusivity(qmeanDh1[1:len(qmeanDh1)-1],dqdt1[0:len(dqdt1)], delt)
print(esth1.params)

estv1=dif.FitDiffusivity(qmeanDz1[1:len(qmeanDz1)-1],dqdt1[0:len(dqdt1)], delt)
print(estv1.params)

print esth1.summary()
print estv1.summary()
print qmeanDh1<jupyter_output>[[  1.10370376e-05]
 [             nan]
 [             nan]
 [             nan]
 [             nan]
 [             nan]
 [             nan]
 [             nan]
 [             nan]
 [             nan]
 [             nan]
 [             nan]
 [             nan]
 [             nan]
 [             nan]
 [             nan]
 [             nan]
 [             nan]
 [             nan]
 [             nan]
 [             nan]]<jupyter_text>Tracer 2<jupyter_code>q2mean2 = dif.CalcVariance(nt, nz, ny, nx, DomVolume, hFacC, drF, rAC, tracer2)
print('done')
dqdt2 = dif.CalcTimeDer(q2mean2, nt, delt)
print('done')
qmeanDh2 = dif.CalcAvgHorGrad(filename2, nt,nz,ny,nx,tracer2,DomVolume,hFacC, drF, rAC)
print('done')
qmeanDz2=dif.CalcAvgVerGrad(filename2, nt,nz,ny,nx,tracer2,DomVolume,hFacC, drF, rAC)
print('done')

esth2=dif.FitDiffusivity(qmeanDh2[1:len(qmeanDh2)-1],dqdt2[0:len(dqdt2)], delt)
print(esth2.params)

estv2=dif.FitDiffusivity(qmeanDz2[1:len(qmeanDz2)-1],dqdt2[0:len(dqdt2)], delt)
print(estv2.params)

print esth2.summary()
print estv2.summary()<jupyter_output><empty_output><jupyter_text>Tracer 3<jupyter_code>q2mean3 = dif.CalcVariance(nt, nz, ny, nx, DomVolume, hFacC, drF, rAC, tracer3)
print('done')
dqdt3 = dif.CalcTimeDer(q2mean3, nt, delt)
print('done')
qmeanDh3 = dif.CalcAvgHorGrad(filename2, nt,nz,ny,nx,tracer3,DomVolume,hFacC, drF, rAC)
print('done')
qmeanDz3=dif.CalcAvgVerGrad(filename2, nt,nz,ny,nx,tracer3,DomVolume,hFacC, drF, rAC)
print('done')

esth3=dif.FitDiffusivity(-qmeanDh3[1:len(qmeanDh3)-1],dqdt3[0:len(dqdt3)], delt)
print(esth3.params)

estv3=dif.FitDiffusivity(-qmeanDz3[1:len(qmeanDz3)-1],dqdt3[0:len(dqdt3)], delt)
print(estv3.params)

print esth3.summary()
print estv3.summary()<jupyter_output><empty_output><jupyter_text>Tracer variance plots:--------------------<jupyter_code>kh1=esth1.params
kv1=estv1.params

kh2=esth2.params
kv2=estv2.params

kh3=esth3.params
kv3=estv3.params



fig1=plt.figure(figsize=(21,3))

plt.subplot(1,3,1)
ax=plt.gca()


plt.plot(qmeanDh1[1:len(qmeanDh1)-1],dqdt1)
plt.plot(qmeanDh2[1:len(qmeanDh2)-1],dqdt2)
plt.plot(qmeanDh3[1:len(qmeanDh3)-1],dqdt3)



plt.xlabel('$<|\partial_{h}q|^2>$')

plt.ylabel('$\partial{<q^2>}/\partial{t}$')
    

plt.subplot(1,3,2)
ax=plt.gca()

plt.plot(qmeanDz1[1:len(qmeanDz1)-1],dqdt1)
plt.plot(qmeanDz2[1:len(qmeanDz2)-1],dqdt2)
plt.plot(qmeanDz3[1:len(qmeanDz3)-1],dqdt3)



plt.xlabel('$<\partial{q^2}/\partial{z}>$')

plt.ylabel('$\partial{<q^2>}/\partial{t}$')

plt.subplot(1,3,3)
ax=plt.gca()


plt.plot(np.arange(0,nt,1),q2mean1)
plt.plot(np.arange(0,nt,1),q2mean2)
plt.plot(np.arange(0,nt,1),q2mean3)




plt.xlabel('hours')

plt.ylabel('$<q^2>$')

plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Numerical diffusivities as a function of prescribed vertical diffusivity<jupyter_code>fig1=plt.figure(figsize=(15,6))

plt.subplot(1,2,1)
ax1=plt.gca()

ax1.plot()
ax1.plot()
#ax.plot((10E-7,10E-3,10E-2,10E-1,1.0,10.0),(esth1.params,esth2.params,esth3.params,esth4.params,esth5.params,esth6.params))
ax1.plot((10**-7,10**-5,10**-4),(kh1,kh2,kh3),'ro')
#ax1.set_yscale('symlog')
#ax1.set_xscale('symlog')
plt.xlabel('Kv')

plt.ylabel('Khnum')
    
plt.subplot(1,2,2)
ax2=plt.gca()

ax2.plot()
ax2.plot()
#ax.plot((10E-7,10E-3,10E-2,10E-1,1.0,10.0),(estv1.params,estv2.params,estv3.params,estv4.params,estv5.params,estv6.params))
ax2.plot((10**-7,10**-5,10**-4),(kv1,kv2,kv3),'bo')

#ax2.set_yscale('symlog')
#ax2.set_xscale('symlog')

plt.xlabel('Kv')

plt.ylabel('Kvnum')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1509303928325293, 'fraction_numerical': 0.05387089363657248, 'mean_word_length': 5.021438450899032, 'pattern_counts': {'":': 0, '<': 50, '<?xml version=': 0, '>': 50, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '530250', 'license': 'Apache-2.0', 'n_tokens_mistral': 3861, 'n_tokens_neox': 3593, 'n_words': 620, 'path': 'NumDiffSensitivityTestGMREDIHydRun01.ipynb', 'repo_name': 'UBC-MOAD/outputanalysisnotebooks', 'size': '8706'}