[ecco-support] question about calc_barostream.m

[Martha Buckley]

I Gael,
Thanks for the reply---I was being silly forgetting about there being an
arbitrary constant.

I did try your new function, and now I get the same result if I integrate
the layer streamfunction as the barotropic streamfunction calculated from
the 3D field.  It is nice that the method of determining the constant now
yields the same result, even if the constant is arbitrary.  Thanks for
modifying the function!

I actually hadn't even noticed that you had sent the other e-mails to
ecco-support, as the reply was just to you.  Thanks for fixing the list
setup.
Martha


On Thu, Mar 17, 2016 at 4:01 PM, gael forget <gforget at mit.edu> wrote:

> 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
>
>
>
> _______________________________________________
> 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/20160318/43ea4a63/attachment-0001.html