[MITgcm-devel] [MITgcm-support] numerical issues with exf_interp_uv
Longjiang Mu
mulongjiang at qq.com
Thu May 25 06:49:50 EDT 2017
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
-------------- next part --------------
A non-text attachment was scrubbed...
Name: uwind_after_interp.png
Type: image/png
Size: 11243 bytes
Desc: not available
URL: <http://mitgcm.org/pipermail/mitgcm-devel/attachments/20170525/60142a34/attachment-0001.png>
More information about the MITgcm-devel
mailing list