Open MaceKuailv opened 9 months ago
BTW, it does not make a difference if one includes hFacS and hFacW
ds['utrans'] = (ds['UVELSTAR']*ds['drF']*ds['dyG']*ds['hFacW']).fillna(0)
ds['vtrans'] = (ds['VVELSTAR']*ds['drF']*ds['dxG']*ds['hFacS']).fillna(0)
I tried to use xgcm for the same calculation,
strm = xr.open_dataset('OCEAN_BOLUS_STREAMFUNCTION_day_mean_1992-01-01_ECCO_V4r4_native_llc0090.nc')
strmx = strm.GM_PsiX.fillna(0)
strmy = strm.GM_PsiY.fillna(0)
u = xgcmgrd.diff(strmx,"k", boundary = 'fill', fill_value = 0.0)/ds.drF
v = xgcmgrd.diff(strmy,"k", boundary = 'fill', fill_value = 0.0)/ds.drF
vstrmx = strmx*ds.dyG
vstrmy = strmy*ds.dxG
xy_diff = xgcmgrd.diff_2d_vector(
{"i": vstrmx, "j": vstrmy}, boundary="fill", fill_value=0.0
)
x_diff = xy_diff["i"]
y_diff = xy_diff["j"]
hDiv = (x_diff + y_diff)
w = hDiv/ds.rA
When I plug this into the divergence calculation above, the divergence is very small. The horizontal components are the same but the vertical velocity at the second level are different.
diffu = (u - ds.UVELSTAR).fillna(0)
diffv = (v - ds.VVELSTAR).fillna(0)
diffw = (w - ds.WVELSTAR).fillna(0)
plt.imshow(diffw[0,1,10].T)
plt.colorbar()
which yields:
@ifenty @owang1 can you help with this? It relates to python toolboxes and ECCOv4r4 which you know better than me
Thank you! Let me know if there is anything that I could help.
Thanks @gaelforget
@MaceKuailv Thank you for the detailed description and the Python code regarding the issue. I will look into it and get back to you.
Thank you! Let me know if there is anything that I could help.
@MaceKuailv sorry that you didn't get a reply sooner. From what I understand, the issue likely comes from either how the vertical bolus velocity is calculated or how it is prepared/distributed.
You mentioned that the issue only appears ... " at grids that have exactly 2 layers in depth around land boundaries, and it only appears in the first two layers. " which is a very important clue.
@owang01 does the model calculate these bolus terms or do we calculate it offline?
I think MITgcm does not have uvelstar as an output variable, The ones that are available are GMpsiX and GMpsiY. Here is a python script of mine to calculate the bolus velocity with xgcm. This one does not create a divergent velocity field.
The code is translated from the code in this repo, so I think the MATLAB tool box probably does not have the same bug.
I think MITgcm does not have uvelstar as an output variable, The ones that are available are GMpsiX and GMpsiY. Here is a python script of mine to calculate the bolus velocity with xgcm. This one does not create a divergent velocity field.
So one can calculate divergence-free bolus velocities from the GMxx output fields from the model, but it seems the code we use has a bug that appears in the case of exactly two depth levels next to land. Is that your understanding?
Yes, exactly.
We calculated these bolus terms offline. It is likely that there was a bug in calculating these bolus terms. I was able to reproduce the issue that @MaceKuailv found but did not find the root cause of the the bug at that time. I will take one more crack at it.
We calculated these bolus terms offline. It is likely that there was a bug in calculating these bolus terms. I was able to reproduce the issue that @MaceKuailv found but did not find the root cause of the the bug at that time. I will take one more crack at it.
@owang01 Mace seems to calculate the w that yields a zero divergence. I don't fully understand the calculation (should hFac* be involved?) because it seems zero horizontal divergence yields a zero w. Shouldn't it be zero horizontal divergence yields zero vertical divergence, wk = w{k+1}?
Hi!
I am using the ECCO v4r4 product and realized that the bolus velocity is not purely divergent-free. I think this is the repo that hosts the code that is responsible for converting GMPsiX and GMPsiY to UVELSTAR, VVELSTAR, WVELSTAR (https://github.com/MITgcm/gcmfaces/blob/master/gcmfaces_calc/calc_bolus.m#L1).
The finite divergence appears at grids that have exactly 2 layers in depth around land boundaries, and it only appears in the first two layers. This is what the convergence of bolus velocity looks like:
Here is my code to go from downloading to producing the above plot. (Excuse me for not able to write in MATLAB)
For comparison, the divergence should be many orders of magnitude smaller. Here is what it looks like in the third layer.
The code for reproducing this is:
Although I don't know how to write MATLAB, please let me know how I can help.