Good day everyone . I have some difficulties to obtain geostrophic current by using ROMS output. I am using the script provided by John Wilkin (uvz_from_clim_stdlev.m) to calculate geostrophic current with the use of climatological ROMS output. Unfortunately, this script seem to only compatible with old version of ROMS which using z coordinate (constant depth for each layer throughout lon and lat) instead of the sigma coordinate (different depth for each layer throughout lon and lat) as in the latest ROMS version. I try to edit several lines of the script (see attachment and indicated by "Edited by Ernest") so it able to work with latest version of ROMS, which including the section of:-
extract pressure:
pres = squeeze(z(:,J,:));
pres = flipud(abs(pres));
v component velocity:
ga0 = repmat(-9.81*zeta',size(z,1),1);
u component velocity:
ga0 = repmat(-9.81*zeta,size(z,1),1);
But, I have no idea on how to edit command under the section of "relative to input pref", it keeps to return me "Index in position 2 exceeds array bounds".
l_ref = min(find(pres >= abs(pref)));
hclm(mon,:,:) = squeeze(dha(mon,l_ref,:,:)-dha(mon,1,:,:));
Any idea of it or any latest version of the script of uvz_from_clim_stdlev.m so that I can use for the latest version of ROMS output?
Kindly find the attachment of the script that I have edited.
My data have 12 time frame and 26 sigma level.
p/s: Your help is much appreciated.
Calculation of Geostrophic Current
Calculation of Geostrophic Current
- Attachments
-
- uvz_from_clim_stdlev.m
- (5.94 KiB) Downloaded 464 times
Re: Calculation of Geostrophic Current
That is a script designed to do what it says in the help:
It does not operate on ROMS output. If it were a function designed specifically for ROMS I would have called it roms_(something).m.
There has never been an "old version of ROMS" using a z coordinate.
You'll have to write the script you want from first principles unless someone out there has a Matlab function that does what you want. I do not.
Code: Select all
% Compute u,v geostrophic from a standard level t,s,ssh climatology
There has never been an "old version of ROMS" using a z coordinate.
You'll have to write the script you want from first principles unless someone out there has a Matlab function that does what you want. I do not.
John Wilkin: DMCS Rutgers University
71 Dudley Rd, New Brunswick, NJ 08901-8521, USA. ph: 609-630-0559 jwilkin@rutgers.edu
71 Dudley Rd, New Brunswick, NJ 08901-8521, USA. ph: 609-630-0559 jwilkin@rutgers.edu
Re: Calculation of Geostrophic Current
Dear all,
actually, I had a similar question last week and I'm quite happy there is some discussion about that.
However, as this is really a basic issue of physical oceanography, I assume there are other posts concerning this topic, which I might have missed.
My comment is more a fundamental one about the procedure how to obtain these currents.
It might not be a direct answer to the first question.
My present understanding is that the method to use is depending on the type of the vertical coordinate one is using:
a) Using the standard hydrographic data with pressure in dbar as the vertical coordinate, one might calculate the geostrophic currents via gradients of geopotential along a constant pressure surface. This is actually what the matlab script is doing, which you mentioned and attached.
b) However, if your vertical coordinate is somehow related to geopotential (z coordinate), then you might estimate the currents from gradients of pressure on levels of constant z.
So, the roms output, which delivers zeta as well as potential temperature (Tp) and salinity (S) as function of s-coordinates (which can be transformed to z-coordinates) might make the approach b) more natural.
So, when also thinking that approach b) might be suitable, then the pressure at some depth level might be obtained as:
p(z) = p_atmo + g*rho_0*zeta - g*rho_0*z + g* int_(z)^(zeta) rho_ano dz,
where rho_0 is some reference density e.g. 1000 kg/m^3 and rho_ano the density anomaly of in situ density to this reference.
Assuming a homogeneous atmospheric pressure, then the active part making the geostrophic currents might be obtained by gradients of the scaled pressure p* = p / (rho_0*g):
p*(z) = zeta + int_(z)^(zeta) (rho_ano/rho_0) dz
p*(z) = zeta + dh(z)
where the last integral might be called dynamic height dh and is obtained by vertical integration of density anomaly.
So, the first term is due to variations in zeta and the second due to density variations.
Now, the question might be how to calculate the second part from profiles of Tp and S obtained from the roms model.
One has to use a proper equation of state (eos), which might be delivered by the CSIRO Division Oceanography seawater tools.
As an alternative, I found some matlab script within my roms tools called roms_eos.m, which finally estimates the density anomaly of in situ density from profiles of Tp and S and using z as some pressure variable in dbar.
So, I was a bit sceptical if this roms_eos is accurate enough, especially when using z as some kind of pressure.
Therefore, I made a short test script assuming analytical profiles of Tp and S and 1) estimating the dynamical height dh from roms_eos or 2) from the CSIRO Division Oceanography seawater tools together with a 'better' pressure obtained from hydrodynamic relation and sw-eos solving the differential equation for pressure with ode-tools.
Please find the matlab script attached. (in case you want to let it run, be sure you have the CSIRO Division Oceanography seawater tools and the roms_eos.m in your Matlab search path).
In the end, from this simple test, the roms_eos seems accurate enough to estimate a good dh, which can be used together with zeta to further estimate the geostrophic currents.
I really hope, my statements and the script are correct!!!
Best, Karsten
actually, I had a similar question last week and I'm quite happy there is some discussion about that.
However, as this is really a basic issue of physical oceanography, I assume there are other posts concerning this topic, which I might have missed.
My comment is more a fundamental one about the procedure how to obtain these currents.
It might not be a direct answer to the first question.
My present understanding is that the method to use is depending on the type of the vertical coordinate one is using:
a) Using the standard hydrographic data with pressure in dbar as the vertical coordinate, one might calculate the geostrophic currents via gradients of geopotential along a constant pressure surface. This is actually what the matlab script is doing, which you mentioned and attached.
b) However, if your vertical coordinate is somehow related to geopotential (z coordinate), then you might estimate the currents from gradients of pressure on levels of constant z.
So, the roms output, which delivers zeta as well as potential temperature (Tp) and salinity (S) as function of s-coordinates (which can be transformed to z-coordinates) might make the approach b) more natural.
So, when also thinking that approach b) might be suitable, then the pressure at some depth level might be obtained as:
p(z) = p_atmo + g*rho_0*zeta - g*rho_0*z + g* int_(z)^(zeta) rho_ano dz,
where rho_0 is some reference density e.g. 1000 kg/m^3 and rho_ano the density anomaly of in situ density to this reference.
Assuming a homogeneous atmospheric pressure, then the active part making the geostrophic currents might be obtained by gradients of the scaled pressure p* = p / (rho_0*g):
p*(z) = zeta + int_(z)^(zeta) (rho_ano/rho_0) dz
p*(z) = zeta + dh(z)
where the last integral might be called dynamic height dh and is obtained by vertical integration of density anomaly.
So, the first term is due to variations in zeta and the second due to density variations.
Now, the question might be how to calculate the second part from profiles of Tp and S obtained from the roms model.
One has to use a proper equation of state (eos), which might be delivered by the CSIRO Division Oceanography seawater tools.
As an alternative, I found some matlab script within my roms tools called roms_eos.m, which finally estimates the density anomaly of in situ density from profiles of Tp and S and using z as some pressure variable in dbar.
So, I was a bit sceptical if this roms_eos is accurate enough, especially when using z as some kind of pressure.
Therefore, I made a short test script assuming analytical profiles of Tp and S and 1) estimating the dynamical height dh from roms_eos or 2) from the CSIRO Division Oceanography seawater tools together with a 'better' pressure obtained from hydrodynamic relation and sw-eos solving the differential equation for pressure with ode-tools.
Please find the matlab script attached. (in case you want to let it run, be sure you have the CSIRO Division Oceanography seawater tools and the roms_eos.m in your Matlab search path).
In the end, from this simple test, the roms_eos seems accurate enough to estimate a good dh, which can be used together with zeta to further estimate the geostrophic currents.
I really hope, my statements and the script are correct!!!
Best, Karsten
- Attachments
-
- test_dynamic_height_roms_eos.m
- (7.63 KiB) Downloaded 519 times
-
- Posts: 4
- Joined: Tue Oct 06, 2020 8:02 pm
- Location: Department of Geophysics Banaras Hindu University
Re: Calculation of Geostrophic Current
Hi Karsten,
Hope you are doing well.
I wish to initiate again the previous geostrophic current issue for ocean ROMS community. I am also try to calculate it in Arabian Gulf with 1/120 degree resolution output. I posted one actual problem regarding this on ROMS discussion recently. I just curious how we can get gradient of pressure on skew curvilnear grid of ROMS output without using any in built ROMS or other function. My intention is to calculate it with mathematical formula. I did try to get it with by my own mathematical calculation but still getting a lot of differences in actual and geostropic current which seems not good as usual. I have both data on sigma as well as on Z co-ordinate but the the issue is only how I can take gradient of pressure field on this non-symmetric lat-lon coordinate.
Anurag
Hope you are doing well.
I wish to initiate again the previous geostrophic current issue for ocean ROMS community. I am also try to calculate it in Arabian Gulf with 1/120 degree resolution output. I posted one actual problem regarding this on ROMS discussion recently. I just curious how we can get gradient of pressure on skew curvilnear grid of ROMS output without using any in built ROMS or other function. My intention is to calculate it with mathematical formula. I did try to get it with by my own mathematical calculation but still getting a lot of differences in actual and geostropic current which seems not good as usual. I have both data on sigma as well as on Z co-ordinate but the the issue is only how I can take gradient of pressure field on this non-symmetric lat-lon coordinate.
Anurag