[MITgcm-devel] [MITgcm-support] numerical issues with exf_interp_uv
Jean-Michel Campin
jmc at mit.edu
Fri May 26 14:02:47 EDT 2017
Hi Longjiang,
Things are not (yet) completely clear to me, for instance
> It's very interesting that the data at 90N is almost the same as
> that at 89.7222N.
but from the previous email, I thought that the input files have
values every 0.55555:
> uwind_lat0 = 30.0006000000D0,
> uwind_lat_inc = 108*0.5555500D0,
which would mean that the 2 Northernmost lines are at 89.44385 N and 89.9994 N
Anyway, in the first attempt, when you specify
> uwind_lat0 = 30.0006000000D0,
> uwind_lat_inc = 108*0.5555500D0,
there is a potential for having large gradient near the pole if the Northernmost
(~ at the N.pole) is not consistent with single wind vector.
But at this point, I would not conclude that exf_interp_uv.F is doing
the wrong thing with bi-cubic interp.
The solution would be:
1) use bi-linear interp (which does not try to preserve gradient continuity).
- or -
2) fix the input data set by replacing the Northernmost line with
432 values that are consistent with a single wind vector, so that
there is no more unrealistic large gradient in the input field.
3) An other option would be to discard this Northernmost line;
but might not be convenient + might run into the other problem below.
Regarding your second attempt with:
> uwind_lat0 = 30.0006000000D0,
> uwind_lat_inc = 107*0.5555500D0,0.27777475D0,
the problem is different because the input field might not contain
large gradient but exf_interp_uv.F introduces artificial large gradient
and this need to be fixed.
You may want to try:
> uwind_lat0 = 30.0006000000D0,
> uwind_lat_inc = 107*0.5555500D0,0.28D0,
which will results in initial y_in(nyIn+1) > 90. that will be changed to 90.N
(like in the original first attempt) and no artificial gradient from
exf_interp_uv.F
Cheers,
Jean-Michel
On Thu, May 25, 2017 at 06:49:50PM +0800, Longjiang Mu wrote:
> Hi Jean-Michel,
>
> Sorry for my late reply, trying to find some time for test.
>
> 1. about the input data
> I have checked over the input data again.
> The data is downloaded from ECMWF with data range in 30 -90N and
> grids as archive.
> It's very interesting that the data at 90N is almost the same as
> that at 89.7222N.
> So, as you mentioned, large gradients can be found around north pole.
> That definitely will cause the blow up at the north pole.
>
> 2. The forcing data at 90N should be the same, or with very small
> deviations.
> But I think it's not the main reason caused the extreme values.
> Please see 3 below.
>
> 3. about the lat, lon defined in data.exf
> As your reminding, I change the *wind_lat_inc to be
> (107*0.5555500D0,0.27777475D0), which were (108*0.5555500D0) before:
> >uwind_lon0 = 0.83261501D0,
> > uwind_lon_inc = 0.833333333299D0,
> > uwind_lat0 = 30.0006000000D0,
> > uwind_lat_inc = 107*0.5555500D0,0.27777475D0,
> > uwind_nlon = 432,
> > uwind_nlat = 109,
> >#uwind_interpMethod = 1,
> >
> > vwind_lon0 = 0.83261501D0,
> > vwind_lon_inc = 0.833333333299D0,
> > vwind_lat0 = 30.0006000000D0,
> > vwind_lat_inc = 107*0.5555500D0,0.27777475D0,
> > vwind_nlon = 432,
> > vwind_nlat = 109,
> >#uwind_interpMethod = 1,
>
> Unfortunately the model still blew up finally.
> This is because in exf_interp_uv.F,
> >128 C-- setup input latitude grid
> >129 y_in(1) = lat_0
> >130 DO j=1,nyIn+1
> >131 i = MIN(j,nyIn-1)
> >132 y_in(j+1) = y_in(j) + lat_inc(i)
> >133 ENDDO
> >134 y_in(0) = y_in(1) - lat_inc(1)
> >135 y_in(-1)= y_in(0) - lat_inc(1)
> , in our case,
> >y_in(nyIn)=y_in(nyIn-1)+0.27777475D0=89.72222475
> >y_in(nyIn+1)=y_in(nyIn)+0.27777475D0=89.99999950
> , the y_in(nyIn+1) is very close to 90.0 but still satisfies
> >j = nyIn+2
> >IF ( ABS(y_in(j-1)).LT.yPole .AND. ABS(y_in(j)).GT.yPole ) THEN
> > y_in(j) = yPole
> Then , y_in becomes:
>
> >y_in(nyIn): 89.72222475
> >y_in(nyIn+1): 89.99999950
> >y_in(nyIn+2): 90.00000000
>
> , and before are:
>
> >y_in(nyIn): 89.999998569488525
> >y_in(nyIn+1): 90.000000000000000
> >y_in(nyIn+2): 90.000001430511475
>
>
> So, if the y_in passed into the LAGRAN function,
>
> > denom = denom*(a(i) - a(k))
> > numer = numer*(x - a(k))
>
> numer/denom can be 10^3 .
> What we learned above is that the lat(end) + lat_inc(end) should
> never be too close to 90.
> The figure attached is the extreme values near north pole (the model
> domain: ECCO2_panArctic_18km(420x384)).
>
> Best Regards,
> Longjiang
>
> ??? 2017???05???24??? 05:50, Jean-Michel Campin ??????:
> >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
>
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel
More information about the MITgcm-devel
mailing list