[MITgcm-devel] [MITgcm-support] numerical issues with exf_interp_uv
Longjiang Mu
mulongjiang at qq.com
Sun May 28 05:44:01 EDT 2017
Hi Jean-Michel,
Thanks for your suggestions.
I have not explained the details about the data.
The data was downloaded from TIGGE in ECMWF with last 3 latitudes
88.6111, 89.1667, 89.7222 as default.
To save space, I just downloaded the data within 30N-90N that the domain
can be selected during retrieving.
The I found the last 3 latitudes was 88.8889, 89.4444, 90.0000.
The ECMWF has done the interpolation to my grids, but the data at 90N(my
version) is the same as 89.7222(the default version).
So..., I changed the last lat_inc. Then the model blew up at the first step.
I also used a slight bigger lat_inc to cheat the model like your
suggestions, it works well.
Best,
Longjiang
在 2017年05月27日 02:02, Jean-Michel Campin 写道:
> 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