[MITgcm-support] Re: doubly-reentrant

[Baylor Fox-Kemper]

Yes Joe (and mitgcm-support for posterity),
   That's true.  As I understand the finite volume part of it, you 
figure the flux through every face of a volume centered on the 
(temperature, salt, divergence) point.

First, you also don't change tracer inside topography:  in 
gad_advection, when you figure the flux divergence you get lines like:

            localTij(i,j)=localTij(i,j)-dTtracerLev(k)*
      &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
      &       *recip_rA(i,j,bi,bj)
      &       *( af(i+1,j)-af(i,j)
      &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
      &        )

The hFacC (from ini_masks_etc.F) here gives you the topography, i.e., 
in a wall you don't change temp.  But, you might still have a flux 
through topography...

There is a mask hfacW that is 0 for blocked western cell faces, 1 for 
open cell faces, and a fraction for lopped cell faces.  Similarly hfacS 
is for southern faces...

You then use this factor to calculate cell face area (which you do 
every timestep in calc_common_factors, called from thermodynamics.F to 
be sure that the free surface is treated properly):

C--   Calculate tracer cell face open areas
       DO j=jMin,jMax
        DO i=iMin,iMax
         xA(i,j) = _dyG(i,j,bi,bj)
      &   *drF(k)*_hFacW(i,j,k,bi,bj)
         yA(i,j) = _dxG(i,j,bi,bj)
      &   *drF(k)*_hFacS(i,j,k,bi,bj)
        ENDDO
       ENDDO

So, now xa has info about the topography built in.

For the advective fluxes, you also make sure topography is built into 
uTrans and vTrans (also from calc_common_factors):

C--   Calculate velocity field "volume transports" through
C--   tracer cell faces.
       DO j=jMin,jMax
        DO i=iMin,iMax
         uTrans(i,j) = uVel(i,j,k,bi,bj)*xA(i,j)
         vTrans(i,j) = vVel(i,j,k,bi,bj)*yA(i,j)
        ENDDO
       ENDDO

So for diffusive fluxes, you find in gad_diff_x that xA comes in to set 
the b.c.:

       DO j=1-Oly,sNy+Oly
        dfx(1-Olx,j)=0.
        DO i=1-Olx+1,sNx+Olx
         dfx(i,j) = -diffKh*xA(i,j)
      &      *_recip_dxC(i,j,bi,bj)
      &      *(Tracer(i,j)-Tracer(i-1,j))
      &      *CosFacU(j,bi,bj)
        ENDDO
       ENDDO

And for all of the various flux limiting advection schemes 
(gad_*_adv_*.F) you find that the fluxes are always proportional in the 
end to a uFlux which already has the b.c. built in...or that there is 
an xA to make it work.

Similar detective work looking for hFacUp will tell you when vertical 
fluxes are blocked by topography...

Cheers,
    -Baylor


On Mar 31, 2005, at 8:51 AM, Joseph LaCasce wrote:

>
> Hey Baylor---
>
> What does that mean, exactly? If the boundary lies, say, in the
> middle of a channel, you still have to impose no flux on either
> side of it, don't you?
>
> j
>
> -- 
> J.H. LaCasce, Phone: +47 22 96 33 11; FAX: +47 22 69 63 55
> Norwegian Meteorological Institute, P.O. Box 43 Blindern, 0313 Oslo
> or (UPS) Niels Henrik Abelsvei 40, 0371 Oslo, Norway
> http://met.no/english/r_and_d_activities/people/josephhl.html
>
>
>