[Thread Prev][Thread Next][Index]
Re: 2D varying z-axis
Hi again,
I solved the remaining problem regarding the partial filled grid cells.
In principle I shifted the calculation of the grid cells from central
points to interfaces.
The new script is:
use test
! calculate delta z and interface depths
LET dz=zbox[g=temp]
let interface=dz[k=@rsum]
! create mask of complete filled cells
LET zmask = IF interface LT bottom THEN 1
! identify partial filled cell and calculate the ratio=filled/complete
LET margin_mask = IF interface[k=@shf:-1] LT bottom AND interface GT
bottom THEN (bottom-interface[k=@shf:-1])/dz
! merge both masks
LET weighted_mask = MISSING(zmask, margin_mask)
! weight variable with final mask
LET summand = temp * weighted_mask * ZBOX[g=temp]
! integrate
LET result = summand[k=@sum]
By using the following two lines you can test the results with a given
depth and play around with values of bottom. Please make sure that you
are only use values that are IN the partial filled interface (not at the
boundary and not outside), i.e. the deepest grid cell that contains the
variable:
let dum=1+temp*0
list/l=1 dum[z=@din],dum[k=1:26@din]+dum[k=27]*270,result
I hope that's it.
Arne
Steve Hankin wrote:
>
> Hi Arne,
>
> hmmm ... I think there is a flaw in the logic I suggested. The test
> "IF zcoord LT bottom"
> tests the location of the grid point, which is internal to the grid cell. Really the test
> for whether the bottom lies beyond the grid cell should be based on the bathymetry being
> deeper than the BOTTOM of the grid cell. And the computation of the weight for the grid
> cells that are "on the margin" should be based on the TOP edge of the grid cell.
>
> If your Z axis is regular this is a very straightforward correction. If not there are some
> details locating the grid upper/lower boundard ( (zcoord[shifted] + zcoord)/2 is correct
> unless you have explicitly included cell boundaries in your netCDF data set.). I have some
> deadlines I have to meet just now ...
>
> Also try listing the pieces of the calculation in z
> list/k=1:27 temp, weighted_mask, ZBOX[g=temp], summand
> The problem should be apparent from this, I'd think.
>
> In fiddling you might try replacing ZBOX by (zcoord[shifted] - zcoord) as an experiment,
> though this looks OK.
>
> - steve
>
> ======================================
>
> Arne Biastoch wrote:
>
> > Steve Hankin wrote:
> >
> > Steve,
> >
> > > I'd encourage you to run this simple test case ...
> > > e.g.
> > > let myvar1 = 1 + 0 * a_var_on_my_grid ! value is 1.0 everywhere
> > >
> > > and confirm that integrating myvar1 gives you the values of the 2D bathymetry, itself,
> > > to adequate accuracy.
> >
> > good idea. I ran:
> >
> > yes? let dum=1+0*temp
> > yes? LET summand = dum * weighted_mask * ZBOX[g=temp]
> > yes? yes? LET result = summand[k=@sum]
> > yes? yes? list/l=1 dum[z=@din],dum[k=1:26@din]+dum[k=27]*270,result
> > *** NOTE: Ambiguous coordinates on T axis: DUM * WEIGHTED_MASK *
> > ZBOX[G=TEMP]
> > DATA SET: ./test.cdf
> > LONGITUDE: 50W(-50)
> > LATITUDE: 56N
> > TIME: 30-APR-1992 12:00
> > Column 1: DUM[Z=0:6525@DIN] is 1+0*TEMP
> > Column 2: EX#2 is DUM[K=1:26@DIN]+DUM[K=27]*270
> > Column 3: RESULT[Z=0:6525] is SUMMAND[K=@SUM]
> > DUM EX#2 RESULT
> > I / *: 4100. 3745. 3436.
> >
> > But the bottom is at 3745. So my calculation by hand (value 2 is
> > correct), but the one that is calculated using the weighted mask is even
> > lower than down to the level one point above the bottom:
> > yes? yes? list/l=1 dum[k=1:26@din]
> > 1+0*TEMP
> > LONGITUDE: 50W(-50)
> > LATITUDE: 56N
> > DEPTH (m): 0 to 3475 (integrated)
> > TIME: 30-APR-1992 12:00
> > DATA SET: ./test.cdf
> > 3475.
> > yes?
> >
> > Any suggestions?
> >
> > Ciao,
> > Arne
> >
> > > =======================
> > >
> > > Arne Biastoch wrote:
> > >
> > > > Hi Ferreters,
> > > >
> > > > Steve (Hankin) gave me a solution to the problem which worked right from
> > > > the start (I expected this when I got the answer from the master
> > > > itself...). But also thanks to Graham Gladman who gave me a lot of
> > > > suggestions.
> > > >
> > > > Before I come to the solution I should mention that the value is
> > > > slightly different from what I calculated by hand (compare values 2 and
> > > > 3 in the list command below). I currently don't have a clue why. Maybe
> > > > it has to do with details how FERRET integrates.
> > > >
> > > > Here are Steves lines:
> > > > The simplistic answer is to use a vertical mask:
> > > > LET zcoord = Z[g=myvar]
> > > > LET zmask = IF zcoord LT bathymetry2Dfield THEN 1
> > > > Then apply this mask to the quantity in question before integration. The
> > > > effect of the approach is to count every data point either 100% or 0%
> > > >
> > > > It sounds like you want to weight the bottom points by where within the
> > > > grid cell the bottom lies. For that you could use a "weighted mask"
> > > > LET margin_mask = IF zcoord[k=@shf] GT bathymetry2Dfield AND zcoord
> > > > GT
> > > > bathymetry2Dfield THEN
> > > > (bathymetry2Dfield-zcoord)/(zcoord[k=@shf]-zcoord)
> > > > LET weighted_mask = MISSING(zmask, margin_mask)
> > > >
> > > > Then
> > > > LET summand = myvar * weighted_mask * ZBOX[g=myvar]
> > > > LET result = summand[k=@sum]
> > > >
> > > > I ran this on a testfile (attached to this email) and got:
> > > >
> > > > NOAA/PMEL TMAP
> > > > Program FERRET
> > > > Version 5.22 - 07/27/00
> > > > 27-Sep-00 10:44
> > > >
> > > > yes? use test
> > > > yes? show data
> > > > currently SET data sets:
> > > > 1> ./test.cdf (default)
> > > > name title I J
> > > > K L
> > > > TEMP POTENTIAL TEMPERATURE 1:1 1:1 1:30
> > > > 1:1
> > > > BOTTOM (-1)*DEPTH[D=depth] 1:1 1:1 ...
> > > > 1:1
> > > >
> > > > yes? LET zcoord = Z[g=temp]
> > > > yes? LET zmask = IF zcoord LT bottom THEN 1
> > > > yes? LET margin_mask = IF zcoord[k=@shf] GT bottom AND zcoord GT bottom
> > > > THEN (bottom-zcoord)/(zcoord[k=@shf]-zcoord)
> > > > yes? LET weighted_mask = MISSING(zmask, margin_mask)
> > > > yes? LET summand = temp * weighted_mask * ZBOX[g=temp]
> > > > yes? LET result = summand[k=@sum]
> > > > yes?
> > > > yes? list/l=1 temp[z=@din],temp[k=1:26@din]+temp[k=27]*270,result
> > > > *** NOTE: Ambiguous coordinates on T axis: TEMP * WEIGHTED_MASK *
> > > > ZBOX[G=TEMP]
> > > > DATA SET: ./test.cdf
> > > > LONGITUDE: 50W(-50)
> > > > LATITUDE: 56N
> > > > TIME: 30-APR-1992 12:00
> > > > Column 1: TEMP[Z=0:6525@DIN] is POTENTIAL TEMPERATURE (deg C)
> > > > Column 2: EX#2 is TEMP[K=1:26@DIN]+TEMP[K=27]*270
> > > > Column 3: RESULT[Z=0:6525] is SUMMAND[K=@SUM]
> > > > TEMP EX#2 RESULT
> > > > I / *: 13545. 12872. 12286.
> > > > yes?
> > > >
> > > > Ciao,
> > > > Arne
> > > >
> > > > Arne Biastoch wrote:
> > > > >
> > > > > Hi Ferreters,
> > > > >
> > > > > sorry, I have to bother you again. Steves suggestion did not solved by
> > > > > case. But I am sure that this is a general problem that all modelers
> > > > > with variable bottom and free surface formulations should also have.
> > > > >
> > > > > I hope that there might be solutions around (Steve might already have 3
> > > > > lines in mind that do this...).
> > > > >
> > > > > Here's the case:
> > > > >
> > > > > I have a standard netCDF file with one depth axis. Here I would like to
> > > > > integrate not over the full water column but only down to a certain
> > > > > depth given by a 2D field (the bottom is not at the edges of the grid
> > > > > cells).
> > > > >
> > > > > standard case shallower deeper
> > > > > +---+ +---+ +---+
> > > > > | | | | | |
> > > > > | T | | T | | T |
> > > > > | | | | | |
> > > > > +---+ +---+ +---+
> > > > > | | | | | |
> > > > > | T | | T | | T |
> > > > > | | | | | |
> > > > > +---+ +---+ +---+
> > > > > | | | | | |
> > > > > | | | | | |
> > > > > | T | | T | | T |
> > > > > | | | | | |
> > > > > | | |...| <--- | |
> > > > > +---+ <--- +---+ +---+
> > > > > |...| |...| | |
> > > > > |...| |...| |...| <---
> > > > > |...| |...| | T |
> > > > > |...| |...| |...|
> > > > > |...| |...| |...|
> > > > > +---+ +---+ +---+
> > > > >
> > > > > T are the values that are NOT missing values, "<--" marks the bottom,
> > > > > given by a 2D field
> > > > >
> > > > > The most intuitive solution would be:
> > > > >
> > > > > shade T[z=0:depth] whereby depth[i,j]
> > > > >
> > > > > but as I understand FERRET it is not possible to use variable
> > > > > integration limits.
> > > > >
> > > > > I appreciate any help.
> > > > >
> > > > > Ciao,
> > > > > Arne
> > > > >
> > > > > Steve Hankin wrote:
> > > > > >
> > > > > > Hi Arne,
> > > > > >
> > > > > > A good starting point might be the FAQ
> > > > > >
> > > > > > "How do I handle sigma coordinate output in Ferret?"
> > > > > >
> > > > > > at
> > > > > > http://www.ferret.noaa.gov/Ferret/FAQ/data_management/sigma_coordinate_demo.html
> > > > > >
> > > > > > Happy Ferreting - steve
> > > > > >
> > > > > --
> > > > >
> > > > > Dr. Arne Biastoch
> > > > > Scripps Institution of Oceanography phone: +1-858-822-3787
> > > > > Physical Oceanography Research Division fax : +1-858-534-9820
> > > > > MS: 0230 email: abiastoch@ucsd.edu
> > > > > 8605 La Jolla Shores Dr.
> > > > > La Jolla, CA 92093-0230, U.S.A. http://www.ecco.ucsd.edu/~biastoch
> > > >
> > > > --
> > > >
> > > > Dr. Arne Biastoch
> > > > Scripps Institution of Oceanography phone: +1-858-822-3787
> > > > Physical Oceanography Research Division fax : +1-858-534-9820
> > > > MS: 0230 email: abiastoch@ucsd.edu
> > > > 8605 La Jolla Shores Dr.
> > > > La Jolla, CA 92093-0230, U.S.A. http://www.ecco.ucsd.edu/~biastoch
> > > >
> > > > ------------------------------------------------------------------------
> > > > Name: test.cdf
> > > > test.cdf Type: Channel File (application/x-cdf)
> > > > Encoding: base64
> > >
> > > --
> > > Steve Hankin
> > > NOAA/PMEL, 7600 Sand Point Way NE, Seattle, WA 98115-0070
> > > ph. (206) 526-6080 -- FAX (206) 526-6744
> >
> > --
> >
> > Dr. Arne Biastoch
> > Scripps Institution of Oceanography phone: +1-858-822-3787
> > Physical Oceanography Research Division fax : +1-858-534-9820
> > MS: 0230 email: abiastoch@ucsd.edu
> > 8605 La Jolla Shores Dr.
> > La Jolla, CA 92093-0230, U.S.A. http://www.ecco.ucsd.edu/~biastoch
>
> --
>
> | NOAA/PMEL | ph. (206) 526-6080
> Steve Hankin | 7600 Sand Point Way NE | FAX (206) 526-6744
> | Seattle, WA 98115-0070 | hankin@pmel.noaa.gov
--
Dr. Arne Biastoch
Scripps Institution of Oceanography phone: +1-858-822-3787
Physical Oceanography Research Division fax : +1-858-534-9820
MS: 0230 email: abiastoch@ucsd.edu
8605 La Jolla Shores Dr.
La Jolla, CA 92093-0230, U.S.A. http://www.ecco.ucsd.edu/~biastoch
[Thread Prev][Thread Next][Index]
Dept of Commerce /
NOAA /
OAR /
PMEL /
TMAP
Contact Us | Privacy Policy | Disclaimer | Accessibility Statement