[MITgcm-devel] [MITgcm-support] numerical issues with exf_interp_uv
Jean-Michel Campin
jmc at mit.edu
Tue May 23 17:50:52 EDT 2017
Hi Longjiang,
I am switching to mitgcm-devel to sort out the issue.
1) May be I missed it, but I still don't know where you are interpolated to,
i.e., my first question:
> 1) which type of model grid are you using and where does the "bad" interpolated
> values are found (just at the N.pole ?) ?
2) Since the last latitude is very close to the North pole (~ 90 - 1.e-6),
I would assume (and exf interpolation does) that all the 432 values (in lon
direction) represent the same (or nearly the same) wind vector (i.e., the wind
at the North pole) but would be oriented along the local direction (local eastward
and local northward directions) to get the 432 values:
uwind(1:432,nlat) = U_pole * cos( lon*rad ) + V_pole * sin( lon*rad )
vwind(1:432,nlat) = -U_pole * sin( lon*rad ) + V_pole * cos( lon*rad )
And these 432 value should be very close to the single averaged value (after local
rotation) that pkg/exf computes at the north pole.
If it's not the case, then it means that the wind input field that you provide
contains very large gradient near the north pole and, as a consequence, the
bi-cubic interpolation returns large extreme values.
I would be curious to know if you are in this case.
When I try to estimate how the function LAGRAN works, I don't see a clear problem
if there is no gradient (unique wind vector) near N.pole (the 1.e-6 factor
cancel out).
Cheers,
Jean-Michel
On Tue, May 23, 2017 at 09:35:21PM +0800, Longjiang Mu wrote:
> Hi Jean-Michel,
>
> 1) Parts of data.exf:
>
> uwind_lon0 = 0.83261501D0,
> uwind_lon_inc = 0.833333333299D0,
> uwind_lat0 = 30.0006000000D0,
> #uwind_lat0 = 29.99000000D0,
> uwind_lat_inc = 108*0.5555500D0,
> uwind_nlon = 432,
> uwind_nlat = 109,
>
> vwind_lon0 = 0.83261501D0,
> vwind_lon_inc = 0.833333333299D0,
> #vwind_lat0 = 29.99000000D0,
> vwind_lat0 = 30.0006000000D0,
> vwind_lat_inc = 108*0.5555500D0,
> vwind_nlon = 432,
> vwind_nlat = 109,
>
> If the data.exf is configured as above, then y_in calculated in
> exf_interp_uv.F will firstly result:
>
> y_in(nyIn-1): 89.444448590278625
> y_in(nyIn): 89.999998569488525
> y_in(nyIn+1): 90.555548548698425
> y_in(nyIn+2): 91.111098527908325
>
> , then after these codes in exf_interp_uv.F:
>
> C- Add 2 row @ northern end; if one is beyond N.pole, put one @ N.pole
> j = nyIn+1
> IF ( ABS(y_in(j-1)).LT.yPole .AND. ABS(y_in(j)).GT.yPole ) THEN
> y_in(j) = yPole
> y_in(j+1) = 2.*yPole - y_in(j-1)
> #ifdef ALLOW_DEBUG
> prtPole(3) = 1.
> prtPole(3+1) = 2.
> #endif
> ENDIF
> j = nyIn+2
> IF ( ABS(y_in(j-1)).LT.yPole .AND. ABS(y_in(j)).GT.yPole ) THEN
> y_in(j) = yPole
> #ifdef ALLOW_DEBUG
> prtPole(4) = 1.
> #endif
> ENDIF
>
> y_in will be:
>
> y_in(nyIn-1): 89.444448590278625
> y_in(nyIn): 89.999998569488525
> y_in(nyIn+1): 90.000000000000000
> y_in(nyIn+2): 90.000001430511475
>
> , where y_in(nyIn+1) is set to 90.0 _d 0.
> So that, if the y_in above is used in the EXF_INTERPOLATE.F where
> function LAGRAN is called,
> y_in(nyIn) - y_in(nyIn+1) (which is a(i)-a(k) in function
> LAGRAN)
> is very small, then very large numer/denom will be catched.
> And it is only found at north pole as the south pole has not been
> tested yet.
> But we believe that can also happen in the south pole, because the
> same procedure is applied in the routine(line:141-160, v64a).
>
> 2).
> STDOUT.0000 generated using the data.exf in 1) :
>
> 9574 (PID.TID 0000.0001) EXF_INTERP_READ: opening file: u_2010
> 9575 (PID.TID 0000.0001) EXF_INTERP_READ: opening file: v_2010
> 9576 (PID.TID 0000.0001) EXF_INTERP_UV: fileU="u_2010", rec= 1 ,
> x-Per,P.Sym= T T
> 9577 (PID.TID 0000.0001) S.edge (j=-1,0,1) : proc= 0.0 0.0 0.0, yIn=
> 28.889500 29.445050 30.000600
> 9578 (PID.TID 0000.0001) N.edge (j=+0,+1,+2) proc= 0.0 1.1 2.2, yIn=
> 90.000000 90.000000 90.000000
>
> 9580 (PID.TID 0000.0001) EXF_INTERP_READ: opening file: u_2010
> 9581 (PID.TID 0000.0001) EXF_INTERP_READ: opening file: v_2010
> 9582 (PID.TID 0000.0001) EXF_INTERP_UV: fileU="u_2010", rec= 2 ,
> x-Per,P.Sym= T T
> 9583 (PID.TID 0000.0001) S.edge (j=-1,0,1) : proc= 0.0 0.0 0.0, yIn=
> 28.889500 29.445050 30.000600
> 9584 (PID.TID 0000.0001) N.edge (j=+0,+1,+2) proc= 0.0 1.1 2.2, yIn=
> 90.000000 90.000000 90.000000
>
> 3). It works as Dimitris suggested.
>
> Best regards,
> Longjiang
>
> ??? 2017???05???23??? 17:41, mitgcm-support-request at mitgcm.org ??????:
> >Message: 2
> >Date: Mon, 22 May 2017 17:57:32 -0400
> >From: Jean-Michel Campin <jmc at mit.edu>
> >To: mitgcm-support at mitgcm.org
> >Subject: Re: [MITgcm-support] numerical issues with exf_interp_uv?
> >Message-ID: <20170522215732.GA42187 at ocean.mit.edu>
> >Content-Type: text/plain; charset=iso-8859-1
> >
> >Hi Martin,
> >
> >Just to help finding exactly where the problem is:
> >1) which type of model grid are you using and where does the "bad"
> >interpolated
> > values are found (just at the N.pole ?) ?
> >2) when you set:
> > > exf_debugLel=3
> > and assuming you did compile with #define ALLOW_DEBUG,
> > what are the reported lines in STDOUT from S/R EXF_INTERP_UV
> > (e.g., "x-Per,P.Sym=" and " S.edge" and " N.edge" lines ) ?
> >3) trying bi-linear interp instead of bi-cubic (Dimitris' suggestion)
> > could be useful too (just to narrow down the issue).
> >
> >Cheers,
> >Jean-Michel
> >
> >On Mon, May 22, 2017 at 10:03:06AM -0700, Dimitris Menemenlis wrote:
> >>Hi Martin, by default vector fields use bicubic interpolation:
> >> ustress_interpMethod = 12
> >> vstress_interpMethod = 22
> >> uwind_interpMethod = 12
> >> vwind_interpMethod = 22
> >>(this is in order to avoid artifacts in derivative, e.g., wind
> >>stress curl)
> >>while scalar field use bilinear interpolation:
> >> hflux_interpMethod = 1
> >> sflux_interpMethod = 1
> >> ???
> >>(in order to avoid overshoots and undershoots).
> >>This might explain why you only see problem in vector interpolation.
> >>Don???t know how to fix though. Let???s wait for JM to chime in.
> >>
> >>Dimitris
> >>
> >>>On May 22, 2017, at 7:16 AM, Martin Losch <Martin.Losch at awi.de> wrote:
> >>>
> >>>Hi all,
> >>>
> >>>Longjiang and I found a potential instability in exf_interp_uv:
> >>>when the last latitude < 90degN of the data grid is very close
> >>>to 90degN, e.g. 89.999??? then the interpolated velocity
> >>>fields can exceed 1e6, because in the function LAGRAN, the
> >>>variable denom becomes very small (a = py_ind, x = yG):
> >>> denom = denom*(a(i) - a(k))
> >>> numer = numer*(x - a(k))
> >>>so that numer/denom becomes very large. The problem can be
> >>>fixed by making sure the that last grid point is far enough
> >>>from 90degN, but maybe we can find a way of issuing a warning
> >>>or even of stopping the model when abs(py_ind-90deg)<single
> >>>precicsion or similar.
> >>>
> >>>This only happens for the vector field but not the scalar
> >>>fields (in our case), we are not quite sure why.
> >>>
> >>>Martin and Longjiang
> >>>
> >>>Martin Losch // Alfred-Wegener-Institut, Helmholtz Zentrum f?r
> >>>Polar- und Meeresforschung
> >>>Postfach 120161, 27515 Bremerhaven, Germany
> >>>mailto:Martin.Losch at awi.de // Tel./Fax: ++49(471)4831-1872/1797
> >>>http://www.awi.de/ueber-uns/organisation/mitarbeiter/martin-losch.html
> >>>
> >>>Please avoid sending me Word or PowerPoint attachments.
> >>>See http://www.gnu.org/philosophy/no-word-attachments.html
> >>>
> >>>
> >>>_______________________________________________
> >>>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
>
>
>
>
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-support
More information about the MITgcm-devel
mailing list