From b2cbb1dbdbaea9abb8731d5bdcc5ccfea05d04d2 Mon Sep 17 00:00:00 2001 From: hajsong Date: Wed, 24 May 2017 16:33:29 -0400 Subject: [PATCH 1/2] adding estNsq.py A python script to estimate N^2 using T and S --- so_box_biogeo/diags/estNsq.py | 75 +++++++++++++++++++++++++++++++++++ 1 file changed, 75 insertions(+) create mode 100644 so_box_biogeo/diags/estNsq.py diff --git a/so_box_biogeo/diags/estNsq.py b/so_box_biogeo/diags/estNsq.py new file mode 100644 index 00000000..2325c9ee --- /dev/null +++ b/so_box_biogeo/diags/estNsq.py @@ -0,0 +1,75 @@ +import numpy as np +from MITgcmutils import rdmds, densjmd95 +from hspython import loadgrid +import matplotlib.pyplot as plt +""" + Estimating Brunt-Vaisala frequency using T and S + In case when N$^2$ is needed but do not have 'DRHODR' saved, + one can still estimate it using T and S +""" +# +# Constants +# +rhoconst = 1.035e3 +g = 9.81 +# +# Load grid files +# (hspython is available at https://github.com/hajsong/hspython) +# +grd = loadgrid('so_box', varname=['XC','YC','RC','hFacC','DRC','RF']) +[nz, ny, nx] = grd.hFacC.shape +# +# Compute the stratification frequency using DRHODR. +# It is defined at the center of the layer. +# +dRHOdr = rdmds('ocestrat', 9, rec=0); # DRHODR is in the first record in "ocestrat" +Nsq = - dRHOdr*g/rhoconst*grd.mskC +# +# Now, estimate Nsq using T and S. +# When computing "drhodr" at the layer interface, density at upper and lower cell +# is computed using the pressure at the interface. +# Nsq is defined at the layer interface. +# +T = rdmds('dynDiag', 9, rec=2) # THETA is in the third record in "dynDiag" +S = rdmds('dynDiag', 9, rec=3) # SALT is in the fourth record in "dynDiag" +Nsq_TS = np.zeros([nz, ny, nx]) +for k in range(1,nz): + press = -rhoconst*g*grd.RF[k]/1e4 # pressure at the interface + # rho at the center of the upper level + urho = densjmd95(S[k-1, :, :],T[k-1, :, :], press) + # rho at the center of the lower level + lrho = densjmd95(S[k, :, :], T[k, :, :], press) + drhodr = (urho-lrho)/grd.DRC[k]*grd.mskC[k, :, :]*grd.mskC[k-1, :, :] + Nsq_TS[k, :, :] = -drhodr*g/rhoconst +# +# Estimating N$^2$ using "RHOAnoma" is not appropriate +# because "RHOAnoma" is computed using the pressure at that level +# +rho = rdmds('ocestrat', 9, rec=1) + rhoconst # RHOAnoma is in the second record +Nsq_ra = np.zeros([nz, ny, nx]) +for k in range(1,nz): + drhodz = (rho[k-1, :, :] - rho[k, :, :])/grd.DRC[k]\ + *grd.mskC[k, :, :]*grd.mskC[k-1, :, :] + Nsq_ra[k, :, :] = -drhodz*g/rhoconst +# +# Check N$^2$ +# +showz = 5 +showx = 5 +scale = 1e5 + +X, Y = np.meshgrid(grd.YC[:, 0], grd.RC[:showz]) + +f, ax = plt.subplots(1, 3, figsize=(16, 4)) + +im = ax[0].contourf(X, Y, Nsq[:showz, :, showx]*scale, np.arange(0, 15.1 , 1), cmap='Reds') +cb = plt.clabel(im,colors='black',fmt='%3.1f') +ax[0].set_title('N$^2$ from DRHODR [x 10$^5$ s$^{-1}$]', color='black', fontsize=15) + +im = ax[1].contourf(X, Y, Nsq_TS[:showz,:,showx]*scale, np.arange(0,15.1,1), cmap='Reds') +cb = plt.clabel(im,colors='black',fmt='%3.1f') +ax[1].set_title('N$^2$ from T and S [x 10$^5$ s$^{-1}$]', color='black', fontsize=15) + +im = ax[2].contourf(X, Y, Nsq_ra[:showz,:,showx]*scale, np.arange(0,15.1,1), cmap='Reds') +cb = plt.clabel(im,colors='black',fmt='%3.1f') +ax[2].set_title('N$^2$, from RHOAnoma [x 10$^5$ s$^{-1}$]', color='black', fontsize=15) From 53748da0980f7d14d62ccbdcb63ceb78d0067fed Mon Sep 17 00:00:00 2001 From: hajsong Date: Wed, 24 May 2017 17:30:49 -0400 Subject: [PATCH 2/2] estNsq update MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit “estNsq.py” is updated and a file “mitgcmgrid.py” is added. --- so_box_biogeo/diags/estNsq.py | 29 +++++++++------- so_box_biogeo/diags/mitgcmgrid.py | 56 +++++++++++++++++++++++++++++++ 2 files changed, 73 insertions(+), 12 deletions(-) create mode 100644 so_box_biogeo/diags/mitgcmgrid.py diff --git a/so_box_biogeo/diags/estNsq.py b/so_box_biogeo/diags/estNsq.py index 2325c9ee..7545f224 100644 --- a/so_box_biogeo/diags/estNsq.py +++ b/so_box_biogeo/diags/estNsq.py @@ -1,11 +1,12 @@ import numpy as np from MITgcmutils import rdmds, densjmd95 -from hspython import loadgrid +from mitgcmgrid import loadgrid import matplotlib.pyplot as plt """ - Estimating Brunt-Vaisala frequency using T and S + Estimating Brunt-Vaisala frequency using T and S. In case when N$^2$ is needed but do not have 'DRHODR' saved, - one can still estimate it using T and S + one can still estimate it using T and S, + and the appropriate equation of state. """ # # Constants @@ -16,19 +17,19 @@ # Load grid files # (hspython is available at https://github.com/hajsong/hspython) # -grd = loadgrid('so_box', varname=['XC','YC','RC','hFacC','DRC','RF']) +grd = loadgrid('../results', varname=['XC','YC','RC','hFacC','DRC','RF']) [nz, ny, nx] = grd.hFacC.shape # # Compute the stratification frequency using DRHODR. -# It is defined at the center of the layer. +# It is defined at the center of the level. # dRHOdr = rdmds('ocestrat', 9, rec=0); # DRHODR is in the first record in "ocestrat" Nsq = - dRHOdr*g/rhoconst*grd.mskC # # Now, estimate Nsq using T and S. -# When computing "drhodr" at the layer interface, density at upper and lower cell -# is computed using the pressure at the interface. -# Nsq is defined at the layer interface. +# When computing "drhodr" at the interface between tracer cells, +# density at upper and lower cell is computed using the pressure at the interface. +# Nsq is defined at the interface. # T = rdmds('dynDiag', 9, rec=2) # THETA is in the third record in "dynDiag" S = rdmds('dynDiag', 9, rec=3) # SALT is in the fourth record in "dynDiag" @@ -58,18 +59,22 @@ showx = 5 scale = 1e5 -X, Y = np.meshgrid(grd.YC[:, 0], grd.RC[:showz]) +Y, Z = np.meshgrid(grd.YC[:, 0], grd.RC[:showz]) f, ax = plt.subplots(1, 3, figsize=(16, 4)) -im = ax[0].contourf(X, Y, Nsq[:showz, :, showx]*scale, np.arange(0, 15.1 , 1), cmap='Reds') +im = ax[0].contourf(Y, Z, Nsq[:showz, :, showx]*scale, np.arange(0, 15.1 , 1), cmap='Reds') cb = plt.clabel(im,colors='black',fmt='%3.1f') ax[0].set_title('N$^2$ from DRHODR [x 10$^5$ s$^{-1}$]', color='black', fontsize=15) +ax[0].set_xlabel('latitude') +ax[0].set_ylabel('depth [m]') -im = ax[1].contourf(X, Y, Nsq_TS[:showz,:,showx]*scale, np.arange(0,15.1,1), cmap='Reds') +im = ax[1].contourf(Y, Z, Nsq_TS[:showz,:,showx]*scale, np.arange(0,15.1,1), cmap='Reds') cb = plt.clabel(im,colors='black',fmt='%3.1f') ax[1].set_title('N$^2$ from T and S [x 10$^5$ s$^{-1}$]', color='black', fontsize=15) +ax[1].set_xlabel('latitude') -im = ax[2].contourf(X, Y, Nsq_ra[:showz,:,showx]*scale, np.arange(0,15.1,1), cmap='Reds') +im = ax[2].contourf(Y, Z, Nsq_ra[:showz,:,showx]*scale, np.arange(0,15.1,1), cmap='Reds') cb = plt.clabel(im,colors='black',fmt='%3.1f') ax[2].set_title('N$^2$, from RHOAnoma [x 10$^5$ s$^{-1}$]', color='black', fontsize=15) +ax[2].set_xlabel('latitude') diff --git a/so_box_biogeo/diags/mitgcmgrid.py b/so_box_biogeo/diags/mitgcmgrid.py new file mode 100644 index 00000000..673edfe2 --- /dev/null +++ b/so_box_biogeo/diags/mitgcmgrid.py @@ -0,0 +1,56 @@ +from MITgcmutils import rdmds +import numpy as np + +def loadgrid(gridpath, region=None, varname=None): + """ + [Function] + grd=loadgrid(gridname, region=None, varname=None): + + [Description] + Read grid files and load them into the 'grd' object + + [Inputs] + gridpath : A path to the grid files. + region : A list defining the boundary for the grid. + [x0,x1,y0,y1] (default is [0,nx,0,ny]) + varname : A list of strings for the grid data to load + + [Output] + grd : An object containing grid data. + + """ + + if varname is None: + varname = ['XC','YC','RAC','DXC','DYC','hFacC','hFacW','hFacS',\ + 'Depth','RC','RF','DRC','DRF','XG','YG','RAZ','DXG','DYG'] + + class grd(object): + for iv, vname in enumerate(varname): + if region is None: + exec('tmpvar = rdmds("'+gridpath+varname[iv]+'")') + tmpvar = tmpvar.squeeze() + exec(varname[iv]+' = tmpvar') + else: + if vname is 'RC' or vname is 'RF' or vname is 'DRC' or vname is 'DRF': + exec('tmpvar = rdmds("'+gridpath+varname[iv]+'")') + else: + exec('tmpvar = rdmds("'+gridpath+varname[iv]+\ + '",region = '+str(region)+')') + tmpvar = tmpvar.squeeze() + exec(varname[iv]+' = tmpvar') + if vname == 'hFacC': + mskC = hFacC.copy() + mskC[mskC==0] = np.nan + mskC[np.isfinite(mskC)] = 1. + if vname == 'hFacW': + mskW = hFacW.copy() + mskW[mskW==0] = np.nan + mskW[np.isfinite(mskW)] = 1. + if vname == 'hFacS': + mskS = hFacS.copy() + mskS[mskS==0] = np.nan + mskS[np.isfinite(mskS)] = 1. + del tmpvar + del grd.iv,grd.vname + + return grd