interpolation error near the ground ( procx,procy>1 and east/.../north_add_in>1 ) – in #10: INT2LM

in #10: INT2LM

<p> I had a look at the source code and all problems seem to be at the same point. </p> <p> inside ‘interp_utilities.f90’ inside the ‘ <span class="caps"> SUBROUTINE </span> interp_q’ </p> <p> around line 740 we find this: </p> DO l2 = kilats, kilate DO l1 = kilons, kilone IF (lframe) THEN IF (.NOT. lmask_lm(l1,l2)) CYCLE ENDIF j1 = i_index(l1,l2) ! first input grid dimension !IF (j1 == -9999) CYCLE j2 = j_index(l1,l2) ! second input grid dimension !Find out which is the nearest input grid point IF (y_wght(l1,l2) &gt; 0.5_ireals) THEN j4 = j2 + 2 ELSE j4 = j2 – 1 <span class="caps"> ENDIF </span> <p> the last IF causes ‘j4’ to be greater than ‘je_in’ at the northpol and smaller than 1 at the southpol. (independent on procx, procy) <br/> the field px doesn´t seem to be bothered and return some value from the memory instead of an error. </p> <p> I changed it to </p> IF (y_wght(l1,l2) &gt; 0.5_ireals) THEN j4 = j2 + 2 IF( j4 &gt; je_in ) THEN j4 = je_in ENDIF ELSE j4 = j2 – 1 IF( j4 &lt; 1 ) THEN j4 = 1 ENDIF ENDIF <p> This is a workaround for my antarctica domain. and “fixes” both(!) errors (at least superficially): the one at the south pole and the other processor-depentend artefact. </p> <p> For the pacific domain that I uploaded/wrote, the vertical strips persist. It seems to be the same problem but in the other direction/dimension. Right below the IF-statment from above in the sourcecode the following values are used: <br/> “j1+2” and “j1-1” </p> pxj1 = aeastnorth(x_wght(l1,l2), px(j1,j2), px(j1+1,j2), px(j1+2,j2)) <p> or </p> pxj1 = awestsouth(x_wght(l1,l2), px(j1,j2), px(j1+1,j2), px(j1-1,j2)) <p> and these are sometimes also less than 1 or greater than ‘ie_in’, when I run int2lm on multiple processors and ‘east_add_in’. </p> <p> I suppose a similar workaround (changing “j1+2” and “j1-1” sometimes) would work here too, but I hope we can find and fix the real error behind it, which seems to be somewhere around splitting the field for many processors and using the “south/…/east_add_in” </p>

  @rolfzentek in #9c75d74

<p> I had a look at the source code and all problems seem to be at the same point. </p> <p> inside ‘interp_utilities.f90’ inside the ‘ <span class="caps"> SUBROUTINE </span> interp_q’ </p> <p> around line 740 we find this: </p> DO l2 = kilats, kilate DO l1 = kilons, kilone IF (lframe) THEN IF (.NOT. lmask_lm(l1,l2)) CYCLE ENDIF j1 = i_index(l1,l2) ! first input grid dimension !IF (j1 == -9999) CYCLE j2 = j_index(l1,l2) ! second input grid dimension !Find out which is the nearest input grid point IF (y_wght(l1,l2) &gt; 0.5_ireals) THEN j4 = j2 + 2 ELSE j4 = j2 – 1 <span class="caps"> ENDIF </span> <p> the last IF causes ‘j4’ to be greater than ‘je_in’ at the northpol and smaller than 1 at the southpol. (independent on procx, procy) <br/> the field px doesn´t seem to be bothered and return some value from the memory instead of an error. </p> <p> I changed it to </p> IF (y_wght(l1,l2) &gt; 0.5_ireals) THEN j4 = j2 + 2 IF( j4 &gt; je_in ) THEN j4 = je_in ENDIF ELSE j4 = j2 – 1 IF( j4 &lt; 1 ) THEN j4 = 1 ENDIF ENDIF <p> This is a workaround for my antarctica domain. and “fixes” both(!) errors (at least superficially): the one at the south pole and the other processor-depentend artefact. </p> <p> For the pacific domain that I uploaded/wrote, the vertical strips persist. It seems to be the same problem but in the other direction/dimension. Right below the IF-statment from above in the sourcecode the following values are used: <br/> “j1+2” and “j1-1” </p> pxj1 = aeastnorth(x_wght(l1,l2), px(j1,j2), px(j1+1,j2), px(j1+2,j2)) <p> or </p> pxj1 = awestsouth(x_wght(l1,l2), px(j1,j2), px(j1+1,j2), px(j1-1,j2)) <p> and these are sometimes also less than 1 or greater than ‘ie_in’, when I run int2lm on multiple processors and ‘east_add_in’. </p> <p> I suppose a similar workaround (changing “j1+2” and “j1-1” sometimes) would work here too, but I hope we can find and fix the real error behind it, which seems to be somewhere around splitting the field for many processors and using the “south/…/east_add_in” </p>
RZ replied

I had a look at the source code and all problems seem to be at the same point.

inside ‘interp_utilities.f90’ inside the ‘ SUBROUTINE interp_q’

around line 740 we find this:

DO l2 = kilats, kilate DO l1 = kilons, kilone IF (lframe) THEN IF (.NOT. lmask_lm(l1,l2)) CYCLE ENDIF j1 = i_index(l1,l2) ! first input grid dimension !IF (j1 == -9999) CYCLE j2 = j_index(l1,l2) ! second input grid dimension !Find out which is the nearest input grid point IF (y_wght(l1,l2) > 0.5_ireals) THEN j4 = j2 + 2 ELSE j4 = j2 – 1 ENDIF

the last IF causes ‘j4’ to be greater than ‘je_in’ at the northpol and smaller than 1 at the southpol. (independent on procx, procy)
the field px doesn´t seem to be bothered and return some value from the memory instead of an error.

I changed it to

IF (y_wght(l1,l2) > 0.5_ireals) THEN j4 = j2 + 2 IF( j4 > je_in ) THEN j4 = je_in ENDIF ELSE j4 = j2 – 1 IF( j4 < 1 ) THEN j4 = 1 ENDIF ENDIF

This is a workaround for my antarctica domain. and “fixes” both(!) errors (at least superficially): the one at the south pole and the other processor-depentend artefact.

For the pacific domain that I uploaded/wrote, the vertical strips persist. It seems to be the same problem but in the other direction/dimension. Right below the IF-statment from above in the sourcecode the following values are used:
“j1+2” and “j1-1”

pxj1 = aeastnorth(x_wght(l1,l2), px(j1,j2), px(j1+1,j2), px(j1+2,j2))

or

pxj1 = awestsouth(x_wght(l1,l2), px(j1,j2), px(j1+1,j2), px(j1-1,j2))

and these are sometimes also less than 1 or greater than ‘ie_in’, when I run int2lm on multiple processors and ‘east_add_in’.

I suppose a similar workaround (changing “j1+2” and “j1-1” sometimes) would work here too, but I hope we can find and fix the real error behind it, which seems to be somewhere around splitting the field for many processors and using the “south/…/east_add_in”

View in channel