[MITgcm-support] calculation of divu

[Jean-Michel Campin]

Hi Camille,

I just want to add one thing to Martin's answer:
If you are using Non-Linear free surface, it's a little bit more complicated:
a) without r* coordiante: the surface level thinkness varies in time
 and need to be accounted for.
b) with r* coordinate, div(u) is not zero, as explained in the r* Ocean Modelling
 paper.

Cheers,
Jean-Michel

On Tue, Apr 12, 2016 at 09:57:05AM +0200, Martin Losch wrote:
> Camille,
> 
> you should get div(u)=0 exactly (in Finite Volume Discretization), because ???w??? is computed from this equation. Also, one of the reasons for using finite volumes in the MITgcm is that the divergence is exact.
> 
> Unless you are using a cartesian grid with constant dx/dy, you will have to use finite volumes to compute the divergence. Finite differences with dxf/dyf will not work. See e.g. model/inc/GRID.h for details for the grid configuration. If you are using partical cells, then you???ll need to add the corresponding factors to drf, e.g. hFacW, hFacS, hFacC. If you use the diagnostics-package, then you can save U/V/WVELMASS, where these factors are already included. This is how you could would compute divergence (without paying attention to the details of array sizes etc) in python:
> 
> dxg = rdmds(???DXG???)
> dyg = rdmds(???DYG???)
> ra = rdmds(???RAC???)
> drf = rdmds(???DRF???)
> hc=rdmds(???hFacC')
> hw=rdmds(???hFacW???)
> hs =rdmds(???hFacS???)
> u=rdmds(???U???)
> v=rdmds(???V???)
> w=rdmds(???W???)
> 
> flwm1=0.
> for k in range(len(drf)-1):
> # horizontal divergence:
> 	flu=u[k,:,:]*hw[k,:,:]*dyg*drf[k]
> 	flv=v[k,:,:]*hs[k,:,:]*dxg*drf[k]
> # this is (flu(j,i+1)-flu(j,i) + flu(j+1,i)-flu(j,i))/ra(j,i)
>         hdiv[k,:,:] = [np.roll(flu,-1,1) - flu + np.roll(flv,-1,0) -flv]/(ra*drf[k]*hc[k,:,:])
> 
> # vertical divergence
> 	flw = w[k,:,:]*ra
> 	vdiv[k,:,:]=[flwm1-flw]/(ra*drf[k]*hc[k,:,:])
> 	flwm1=np.copy(flw)
> 
> 
> or similar. I hope I didn???t make a mistake (I usually do).
> 
> Martin
> 
> 
> > On 11 Apr 2016, at 18:05, Camille Mazoyer <mazoyer at univ-tln.fr> wrote:
> > 
> > Dear all,
> > I would like to calculate div u = 0 in the whole domain and at some zoom levels (subdomain near the coast, where I want div u =0 for neglected terms). I did a matlab script for that with u,v,w, dxf, dyf, drf for input, but i'm not very satisfied of it for several reasons.
> > Is there a nice way to calculate divu with some data from the outputs of the model?
> > 
> > Just to be sure, what value of div u should I expect ? (I except small values, but what is a small value for div u?).
> > 
> > Thank you,
> > Camille Mazoyer
> > 
> > -- 
> > ------------------------------------------
> > Camille Mazoyer
> > Phd Student
> > Mediterranean Institute of Oceanography (MIO)
> > Institut de Mathématiques de Toulon (IMATH)
> > Université de TOULON
> > Bat X - CS 60584
> > 83041 TOULON cedex 9
> > France
> > http://mio.pytheas.univ-amu.fr/
> > http://imath.fr/
> > 
> > 
> > _______________________________________________
> > MITgcm-support mailing list
> > MITgcm-support at mitgcm.org
> > http://mitgcm.org/mailman/listinfo/mitgcm-support
> 
> 
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-support