[MITgcm-support] numerical issues with exf_interp_uv
Longjiang Mu
mulongjiang at qq.com
Tue May 23 09:35:21 EDT 2017
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mitgcm.org/pipermail/mitgcm-support/attachments/20170523/06d7956c/attachment.htm>
More information about the MITgcm-support
mailing list