[ecco-support] question about calc_barostream.m
gael forget
gforget at mit.edu
Thu Mar 17 16:01:34 EDT 2016
Hi Martha,
in each horizontal stream function computation there is one indeterminate constant
due to the fact that the velocity field only depends on the stream function gradients.
So the constant offset between meanPsiBT3 and meanPsiBT is immaterial.
At the end of calc_barostream (now lines 128 to 132) you will find that this constant
is evaluated over land as defined by mygrid.mskC(:,:,1) and subtracted. I just
changed the method used to do this so that constant offsets are subtracted
in a more robust way (I think) by always using the same point over Antarctica.
Works?
Gael
On Mar 17, 2016, at 1:52 PM, Martha Buckley <marthabuckley at gmail.com> wrote:
> Thanks for the clarification.
>
> meanPsiBT3-meanPsiBT is now -3.0631 everywhere. Any ideas?
> thanks,
> Martha
>
>
>
> On Wed, Mar 16, 2016 at 8:13 PM, gael forget <gforget at mit.edu> wrote:
> Hi Martha,
>
> to recover the vertically integrated result you want to change:
>> meanPsiLayer(:,:,iz)=tmp.*mygrid.mskC(:,:,iz);
> to
>> meanPsiLayer(:,:,iz)=tmp;
>
> Works?
>
> I reckon that my explanation from the previous email introduced confusion:
> to recover the vertically integrated result you do NOT "want to mask the output
> of calc_barostream.m at level kk with mygrid.mskC(:,:,kk)” but before plotting the
> stream function at level kk by itself I would make sure to apply "mygrid.mskC(:,:,kk)”
> to avoid displaying land values that don’t mean much. Hopefully this clarified things.
>
> Cheers,
> Gael
>
> On Mar 16, 2016, at 10:40 AM, Martha Buckley <marthabuckley at gmail.com> wrote:
>
>> Hi Gael,
>> I updated the code, and now I don't get the singular matrix warnings anymore. If I integrate u, v and then pass the 2D field to calc_barostream, I get the same result as if I pass the 3-D field to calc_barostream. However, if I pass each layer separately to get "psiLayer", if I integrate psiLayer vertically, I still don't recover the barotropic streamfunction. The result is in general too negative, but there are areas when it is too positive and these areas seem to coincide with shallow areas (coastal areas, mid-Atlantic ridge, ares south of New Zealand).
>>
>> Here is my code:
>> meanU meanV are mean of UVELMASS, VVELMASS
>>
>> load mat/diags_grid_parms.mat
>> load my_mat/meanState.mat meanU meanV
>>
>> meanPsiBT=calc_barostream(meanU, meanV);
>> drf=mk3D(mygrid.DRF,meanU);
>> meanU(isnan(meanU))=0; meanV(isnan(meanV))=0;
>> meanPsiBT2=calc_barostream(sum(meanU.*drf,3), sum(meanV.*drf,3));
>>
>> nz=length(mygrid.RC);
>> meanPsiLayer=zeros(mygrid.XC, nz);
>>
>> for iz=1:nz
>> tmp=calc_barostream(meanU(:,:,iz), meanV(:,:,iz));
>> meanPsiLayer(:,:,iz)=tmp.*mygrid.mskC(:,:,iz);
>> end
>>
>> drf=mk3D(mygrid.DRF, meanPsiLayer);
>> meanPsiBT3=nansum(meanPsiLayer.*drf, 3);
>>
>> ____
>> meanPsiBT and meanPsiBT2 are the same, but meanPsiBT3 is not.
>>
>> thanks,
>> Martha
>>
>>
>> On Fri, Mar 11, 2016 at 9:04 AM, gael forget <gforget at mit.edu> wrote:
>> Hi Martha,
>>
>> should work now. Could you please try again after
>> updating calc_barostream.m and diffsmooth2D_div_inv.m ?
>>
>> As it turns out I should have paid more attention to the singular matrix warning.
>> The issue had to do with the creation of isolated canyons when you go down
>> in the water column, which generates the warning during the inversion of the
>> velocity potential. I think that ideally one may want to solve for a specific velocity
>> potential for each isolated region of the global domain at each level kk — but
>> I did not go that far. The revised code instead solves for the velocity potential
>> over the global domain (that is continuous and bridges over all canyons) defined
>> by the surface ocean mask regardless of kk. So you want to mask the output
>> of calc_barostream.m at level kk with mygrid.mskC(:,:,kk).
>>
>> Cheers,
>> Gael
>>
>>
>> On Mar 8, 2016, at 8:34 PM, Martha Buckley <marthabuckley at gmail.com> wrote:
>>
>>> Hi Gael,
>>> I did try passing the integrated 2-D field to the function and I recover the same as the result of passing the 3-D field. Thanks for the suggestion.
>>>
>>> If I, however, sum the streamfunctions calculated in each layer, I don't get the streamfucntion calculated from the 3-D field. I presume that this is because a different divergent part is removed when dealing with the individual layers.
>>> thanks,
>>> Martha
>>>
>>>
>>> On Wed, Mar 2, 2016 at 10:03 AM, gael forget <gforget at mit.edu> wrote:
>>> Hi Martha,
>>>
>>>> Thanks for the tips. In the 2D case, the dz factors are set to 1 directly in the code, so I guess whether they are applied or not makes little difference. I want my diagnostics in SV/m so they don't depend on the layer thickness, so I just didn't apply any dh factors.
>>>
>>> Correct. I overlooked that the default was dry=1 for 2D cases; so if you do not
>>> specify list_factors then only dxg and dyg are applied and this will give you Sv/m.
>>>
>>>> The results look reasonable to me, although I've never looked at the streamfunction in the deep layers before.
>>>> I did try integrating the streamfunction in each layer to see if I could recover the barotropic streamfunction, calculated by passing the full 3-D field to the same function. The patters look qualitatively the same, but there are differences on the order of 5 Sv in the eastern subpolar gyre where the barotropic streamfunction I get summing the layer-by-layer streamfunctions is too weak (not negative enough).
>>>
>>> Checking that you can do the summation of U/VVELMASS*dz in the vertical yourself
>>> and pass this 2D field (in m^2/s) to calc_barostream.m to recover the reference
>>> result (obtained by passing the 3D U/VVELMASS to calc_barostream.m)
>>> should not be too hard and is probably worth doing.
>>>
>>> Cheers,
>>> Gael
>>>
>>> On Mar 1, 2016, at 10:36 AM, Martha Buckley <marthabuckley at gmail.com> wrote:
>>>
>>>> Thanks Gael. Some follow-up is below.
>>>>
>>>> On Mon, Feb 29, 2016 at 4:19 PM, gael forget <gforget at mit.edu> wrote:
>>>> Hi Martha,
>>>>
>>>> I am transferring this thread to ecco-support since the answers may be of interest to others.
>>>>
>>>>> I've used calc_barostream.m many times to calculate the barotropic streamfuntion. I'd like to also use it to calculate the streamfunction in each layer, of course removing the divergent part of the flow. It seems to me that the function can be used for this "as is", by simply passing u, v in a single depth layer, and then the function will return the streamfunction in Sv/m in that layer.
>>>>>
>>>>> Two questions
>>>>> (1) is this use correct?
>>>>> (2) If I do the above to calculate the streamfunction in each layer, I get a warning from diffsmooth2D_div_inv.m that the matrix is singular to working precision for a number of the depth levels.
>>>>
>>>>
>>>> Yes this should work but you want to make sure that you specify an adequate list of scaling factors:
>>>>
>>>> 1) depending on the output you use dh | dz | hfac may have already been factored in by the MITgcm
>>>> diagnostics package; UVELMASS & VVELMASS for example include hfac whereas UVEL & VVEL do not;
>>>> neither include the horizontal and vertical grid spacing (i.e. they are in m/s) but some other diagnostics may.
>>>> 2) in the 3D case you can omit the list_factors argument of calc_barostream to let it apply the default dh and
>>>> dz factors for you. But in 2D cases you want to apply dz beforehand and you might as well also apply dh
>>>> beforehand too. Once you have converted velocities to m3/s you should be able to call calc_barostream
>>>> with list_factors={} as the fourth argument — in analogy with what's done in diags_set_SEAICE.m.
>>>>
>>>> Thanks for the tips. In the 2D case, the dz factors are set to 1 directly in the code, so I guess whether they are applied or not makes little difference. I want my diagnostics in SV/m so they don't depend on the layer thickness, so I just didn't apply any dh factors.
>>>>
>>>> I would disregard the singular matrix warning unless you see something weird in the result.
>>>>
>>>> The results look reasonable to me, although I've never looked at the streamfunction in the deep layers before.
>>>> I did try integrating the streamfunction in each layer to see if I could recover the barotropic streamfunction, calculated by passing the full 3-D field to the same function. The patters look qualitatively the same, but there are differences on the order of 5 Sv in the eastern subpolar gyre where the barotropic streamfunction I get summing the layer-by-layer streamfunctions is too weak (not negative enough).
>>>>
>>>>
>>>>> I'm trying to look at a section across the subtropical/subpolar gyre boundary, which I want to define as where the barotropic streamfunction is zero. I looked a bit at the gcmfaces programs to see what might be helpful, including gcmfaces_lines_zonal.m and gcmfaces_lines_transp.m I'm still having trouble figuring out the best way to do this. Do you have any suggestions?
>>>>
>>>>
>>>> With regard to your earlier question pasted above (sorry for the delay) I would rather suggest using
>>>> gcmfaces_interp_2d by providing the vectors of longitude and latitude that define the contour
>>>> you like in the barotropic stream function. To identify these points you can do something like:
>>>> mskInt=1*(mygrid.Depth>1000);
>>>> [mskCedge,mskWedge,mskSedge]=gcmfaces_edge_mask(mskInt);
>>>> l=convert2vector(mygrid.XC); L=convert2vector(mygrid.YC); m=convert2vector(mskCedge);
>>>> l=l(find(m)); L=L(find(m));
>>>>
>>>> Two comments:
>>>> - gcmfaces_interp_2d.m is a relatively new capability (or at least I revamped it recently) so
>>>> make sure to use the latest gcmfaces before you try this. On a related note I recently revised
>>>> the gcmfaces examples and documentation (wrt gcmfaces_interp_2d.m see example_interp.m).
>>>> - gcmfaces_edge_mask.m also provides the velocity points (mskWedge,mskSedge) that
>>>> you would need to compute transports across your selected domain edge.
>>>>
>>>> Hope this helps.
>>>>
>>>> Cheers,
>>>> Gael
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> --
>>>> Martha W. Buckley
>>>> marthabuckley at gmail.com
>>>> http://sites.google.com/site/marthabuckley/home
>>>
>>>
>>>
>>>
>>> --
>>> Martha W. Buckley
>>> marthabuckley at gmail.com
>>> http://sites.google.com/site/marthabuckley/home
>>> Hi Martha,
>>>
>>> I am transferring this thread to ecco-support since the answers may be of interest to others.
>>>
>>>> I've used calc_barostream.m many times to calculate the barotropic streamfuntion. I'd like to also use it to calculate the streamfunction in each layer, of course removing the divergent part of the flow. It seems to me that the function can be used for this "as is", by simply passing u, v in a single depth layer, and then the function will return the streamfunction in Sv/m in that layer.
>>>>
>>>> Two questions
>>>> (1) is this use correct?
>>>> (2) If I do the above to calculate the streamfunction in each layer, I get a warning from diffsmooth2D_div_inv.m that the matrix is singular to working precision for a number of the depth levels.
>>>
>>>
>>> Yes this should work but you want to make sure that you specify an adequate list of scaling factors:
>>>
>>> 1) depending on the output you use dh | dz | hfac may have already been factored in by the MITgcm
>>> diagnostics package; UVELMASS & VVELMASS for example include hfac whereas UVEL & VVEL do not;
>>> neither include the horizontal and vertical grid spacing (i.e. they are in m/s) but some other diagnostics may.
>>> 2) in the 3D case you can omit the list_factors argument of calc_barostream to let it apply the default dh and
>>> dz factors for you. But in 2D cases you want to apply dz beforehand and you might as well also apply dh
>>> beforehand too. Once you have converted velocities to m3/s you should be able to call calc_barostream
>>> with list_factors={} as the fourth argument — in analogy with what's done in diags_set_SEAICE.m.
>>>
>>> I would disregard the singular matrix warning unless you see something weird in the result.
>>>
>>>> I'm trying to look at a section across the subtropical/subpolar gyre boundary, which I want to define as where the barotropic streamfunction is zero. I looked a bit at the gcmfaces programs to see what might be helpful, including gcmfaces_lines_zonal.m and gcmfaces_lines_transp.m I'm still having trouble figuring out the best way to do this. Do you have any suggestions?
>>>
>>>
>>> With regard to your earlier question pasted above (sorry for the delay) I would rather suggest using
>>> gcmfaces_interp_2d by providing the vectors of longitude and latitude that define the contour
>>> you like in the barotropic stream function. To identify these points you can do something like:
>>> mskInt=1*(mygrid.Depth>1000);
>>> [mskCedge,mskWedge,mskSedge]=gcmfaces_edge_mask(mskInt);
>>> l=convert2vector(mygrid.XC); L=convert2vector(mygrid.YC); m=convert2vector(mskCedge);
>>> l=l(find(m)); L=L(find(m));
>>>
>>> Two comments:
>>> - gcmfaces_interp_2d.m is a relatively new capability (or at least I revamped it recently) so
>>> make sure to use the latest gcmfaces before you try this. On a related note I recently revised
>>> the gcmfaces examples and documentation (wrt gcmfaces_interp_2d.m see example_interp.m).
>>> - gcmfaces_edge_mask.m also provides the velocity points (mskWedge,mskSedge) that
>>> you would need to compute transports across your selected domain edge.
>>>
>>> Hope this helps.
>>>
>>> Cheers,
>>> Gael
>>>
>>>
>>> _______________________________________________
>>> ecco-support mailing list
>>> ecco-support at mit.edu
>>> http://mailman.mit.edu/mailman/listinfo/ecco-support
>>
>>
>>
>>
>> _______________________________________________
>> ecco-support mailing list
>> ecco-support at mit.edu
>> http://mailman.mit.edu/mailman/listinfo/ecco-support
>>
>>
>>
>>
>> --
>> Martha W. Buckley
>> marthabuckley at gmail.com
>> http://sites.google.com/site/marthabuckley/home
>
>
> _______________________________________________
> ecco-support mailing list
> ecco-support at mit.edu
> http://mailman.mit.edu/mailman/listinfo/ecco-support
>
>
>
>
> --
> Martha W. Buckley
> marthabuckley at gmail.com
> http://sites.google.com/site/marthabuckley/home
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.mit.edu/mailman/private/ecco-support/attachments/20160317/73c36b19/attachment-0001.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: smime.p7s
Type: application/pkcs7-signature
Size: 1843 bytes
Desc: not available
Url : http://mailman.mit.edu/mailman/private/ecco-support/attachments/20160317/73c36b19/attachment-0001.bin
More information about the ecco-support
mailing list