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

in #10: INT2LM

<p> PROBLEM: <br/> ————————- <br/> Int2lm seems to make some interpolation error for the initial field when using n&gt;1 processors. <br/> I use ‘east_add_in = 1, south_add_in = 1’ and this caused the problem. Without both parameters the error doesn´t occure. </p> <p> <span class="caps"> OVERVIEW </span> : <br/> ————————- <br/> I stumbled on this error when I was looking at the variable ‘PS’ from my cclm output. (see problem1.png – the white line is the coast of africa) <br/> Some Pixel values are just wrong. I further recognized, that those pixel changed their location (and the number of pixels) when I changed the Parameter ‘procx’ and ‘procy‘ <br/> I also tracked down that the error accually occured during int2lm. The initialfield (at timestep 0) has this error (see problem2.png). It seems best visible in ‘PP’ and not only near the ground but also in higher levels (see problem3.png – level 20 and 39) </p> <p> [computing int2lm on only 1 processor don´t seem to cause this error ? or may the error just lie at the boundary?] </p> <p> <span class="caps"> SETTINGS </span> : to reproduce <br/> ————————- <br/> I used int2lm form the starterpackage 2015 (but the problem occurs also for an old 1.20 Version I was using before) <br/> The namelist is in the file ‘ <span class="caps"> INPUT </span> ’ (with east_add_in = 1, south_add_in = 1) </p> <p> External Data downloaded from WebPEP: <br/> Model version = <span class="caps"> EXTPAR </span> -2.0.2 <br/> pollon=-160 <br/> pollat=90 <br/> polgam=180 <br/> ie_tot=304 <br/> je_tot=304 <br/> startlon_tot=0 <br/> startlat_tot=-40 <br/> dlon=0.125 <br/> dlat=0.125 <br/> oro=1 <br/> orofilter=0 <br/> landuse=2 <br/> soil=1 <br/> tcl=1 <br/> aot=1 <br/> albedo=2 <br/> ————————- </p>

  @rolfzentek in #59ce6f9

<p> PROBLEM: <br/> ————————- <br/> Int2lm seems to make some interpolation error for the initial field when using n&gt;1 processors. <br/> I use ‘east_add_in = 1, south_add_in = 1’ and this caused the problem. Without both parameters the error doesn´t occure. </p> <p> <span class="caps"> OVERVIEW </span> : <br/> ————————- <br/> I stumbled on this error when I was looking at the variable ‘PS’ from my cclm output. (see problem1.png – the white line is the coast of africa) <br/> Some Pixel values are just wrong. I further recognized, that those pixel changed their location (and the number of pixels) when I changed the Parameter ‘procx’ and ‘procy‘ <br/> I also tracked down that the error accually occured during int2lm. The initialfield (at timestep 0) has this error (see problem2.png). It seems best visible in ‘PP’ and not only near the ground but also in higher levels (see problem3.png – level 20 and 39) </p> <p> [computing int2lm on only 1 processor don´t seem to cause this error ? or may the error just lie at the boundary?] </p> <p> <span class="caps"> SETTINGS </span> : to reproduce <br/> ————————- <br/> I used int2lm form the starterpackage 2015 (but the problem occurs also for an old 1.20 Version I was using before) <br/> The namelist is in the file ‘ <span class="caps"> INPUT </span> ’ (with east_add_in = 1, south_add_in = 1) </p> <p> External Data downloaded from WebPEP: <br/> Model version = <span class="caps"> EXTPAR </span> -2.0.2 <br/> pollon=-160 <br/> pollat=90 <br/> polgam=180 <br/> ie_tot=304 <br/> je_tot=304 <br/> startlon_tot=0 <br/> startlat_tot=-40 <br/> dlon=0.125 <br/> dlat=0.125 <br/> oro=1 <br/> orofilter=0 <br/> landuse=2 <br/> soil=1 <br/> tcl=1 <br/> aot=1 <br/> albedo=2 <br/> ————————- </p>

interpolation error near the ground ( procx,procy>1 and east/.../north_add_in>1 )

PROBLEM:
————————-
Int2lm seems to make some interpolation error for the initial field when using n>1 processors.
I use ‘east_add_in = 1, south_add_in = 1’ and this caused the problem. Without both parameters the error doesn´t occure.

OVERVIEW :
————————-
I stumbled on this error when I was looking at the variable ‘PS’ from my cclm output. (see problem1.png – the white line is the coast of africa)
Some Pixel values are just wrong. I further recognized, that those pixel changed their location (and the number of pixels) when I changed the Parameter ‘procx’ and ‘procy‘
I also tracked down that the error accually occured during int2lm. The initialfield (at timestep 0) has this error (see problem2.png). It seems best visible in ‘PP’ and not only near the ground but also in higher levels (see problem3.png – level 20 and 39)

[computing int2lm on only 1 processor don´t seem to cause this error ? or may the error just lie at the boundary?]

SETTINGS : to reproduce
————————-
I used int2lm form the starterpackage 2015 (but the problem occurs also for an old 1.20 Version I was using before)
The namelist is in the file ‘ INPUT ’ (with east_add_in = 1, south_add_in = 1)

External Data downloaded from WebPEP:
Model version = EXTPAR -2.0.2
pollon=-160
pollat=90
polgam=180
ie_tot=304
je_tot=304
startlon_tot=0
startlat_tot=-40
dlon=0.125
dlat=0.125
oro=1
orofilter=0
landuse=2
soil=1
tcl=1
aot=1
albedo=2
————————-

View in channel
<p> This is an issue that is still not solved in INT2LM. There happen to be spurious unrealistic values when crossing the poles. However, this should only affect the inititialization of <span class="caps"> CCLM </span> . After a certain time the fields should look OK in the <span class="caps"> CCLM </span> output. Do you see these features still in <span class="caps"> CCLM </span> output after several days of simulation? </p>

  @burkhardtrockel in #c54b990

<p> This is an issue that is still not solved in INT2LM. There happen to be spurious unrealistic values when crossing the poles. However, this should only affect the inititialization of <span class="caps"> CCLM </span> . After a certain time the fields should look OK in the <span class="caps"> CCLM </span> output. Do you see these features still in <span class="caps"> CCLM </span> output after several days of simulation? </p>

This is an issue that is still not solved in INT2LM. There happen to be spurious unrealistic values when crossing the poles. However, this should only affect the inititialization of CCLM . After a certain time the fields should look OK in the CCLM output. Do you see these features still in CCLM output after several days of simulation?

<p> I think these are two different issues. The namelist&amp;external parameter I uploaded/wrote are not over the poles but over the pacific (east of -160) and this error still occurs. <br/> Even if I leave out the ‘south_add_in’ and use <span class="caps"> ONLY </span> ‘east_add_in = 1’ the error persists. </p> <p> On the other hand for my domain (antartica) the error can be distinguished by using 3 different ‘procx,procy’ setting. <br/> In the picture (see different_errors.PNG) I marked in black circles the 3 different spots. I think What you meant was the bottom right error (at the southpole) </p> <p> Yes this artifacts get “simulated away”, but plan I simulate only 30h or 36h (6h or 12h spinup) and would like to keep this errors low, so I guess at least for the final runs, I have to run them only on one processor? </p>

  @rolfzentek in #e3ca8c6

<p> I think these are two different issues. The namelist&amp;external parameter I uploaded/wrote are not over the poles but over the pacific (east of -160) and this error still occurs. <br/> Even if I leave out the ‘south_add_in’ and use <span class="caps"> ONLY </span> ‘east_add_in = 1’ the error persists. </p> <p> On the other hand for my domain (antartica) the error can be distinguished by using 3 different ‘procx,procy’ setting. <br/> In the picture (see different_errors.PNG) I marked in black circles the 3 different spots. I think What you meant was the bottom right error (at the southpole) </p> <p> Yes this artifacts get “simulated away”, but plan I simulate only 30h or 36h (6h or 12h spinup) and would like to keep this errors low, so I guess at least for the final runs, I have to run them only on one processor? </p>

I think these are two different issues. The namelist&external parameter I uploaded/wrote are not over the poles but over the pacific (east of -160) and this error still occurs.
Even if I leave out the ‘south_add_in’ and use ONLY ‘east_add_in = 1’ the error persists.

On the other hand for my domain (antartica) the error can be distinguished by using 3 different ‘procx,procy’ setting.
In the picture (see different_errors.PNG) I marked in black circles the 3 different spots. I think What you meant was the bottom right error (at the southpole)

Yes this artifacts get “simulated away”, but plan I simulate only 30h or 36h (6h or 12h spinup) and would like to keep this errors low, so I guess at least for the final runs, I have to run them only on one processor?

<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>

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”

<p> ok; (I think )I found the error. I would be grateful if somebody familiar with the source code could check it. Or should I contact the “Current Code Owner”? </p> <p> Inside the ‘src_decomposition.f90’ inside the ‘ <span class="caps"> SUBROUTINE </span> decompose_cm’, there are local variables ‘nzbounds_e, nzbounds_w, nzbounds_s, nzbounds_n’. Normally these have the value 2, and cause 2 additional pixels on each side of the respective coarse grid (which are needed later on during interpolation). </p> <p> But if ‘east/…/north_add_in’ are used these ‘nzbounds_?’ are set to zero and cause (during the end of the 3rd Section of the subroutine) a too small coarse “sub”-grid. </p> <p> I changed it so that the 2 additional pixels are still taken but afterwards a check is made, that caps the value to (1,ie_in_tot) and (1,je_in_tot) </p> !&lt; ZENTEK nzbounds_e = 2 nzbounds_w = 2 nzbounds_s = 2 nzbounds_n = 2 !ZENTEK&gt; ! looking for southern grid point DO j = 1, je_in_tot IF (latitudes_in(j) &gt;= zminlat) THEN jzstart_in = j-nzbounds_s-1 EXIT ENDIF ENDDO ! looking for northern grid point DO j = 1, je_in_tot IF (latitudes_in(j) &gt; zmaxlat) THEN jzend_in = j+nzbounds_n EXIT ENDIF ENDDO ! looking for western grid point DO i = 1, ie_in_tot IF (longitudes_in(i) &gt;= zminlon) THEN izstart_in = i-nzbounds_e-1 EXIT ENDIF ENDDO ! looking for eastern grid point DO i = 1, ie_in_tot IF (longitudes_in(i) &gt; zmaxlon) THEN izend_in = i+nzbounds_w EXIT ENDIF ENDDO !&lt; ZENTEK IF (jzstart_in &lt; 1) THEN jzstart_in = 1 ENDIF IF (jzend_in &gt; ie_in_tot ) THEN jzend_in = ie_in_tot ENDIF IF (izstart_in &lt; 1) THEN izstart_in = 1 ENDIF IF (izend_in &gt; ie_in_tot) THEN izend_in = ie_in_tot ENDIF !ZENTEK&gt; <p> All in all I think these 4+3*4 lines should solve the problem. </p> <p> Leaving only the small error directly over the south pole, which can be ‘workarounded’(*) by the 3 lines of code from my last post. </p>

  @rolfzentek in #14467c7

<p> ok; (I think )I found the error. I would be grateful if somebody familiar with the source code could check it. Or should I contact the “Current Code Owner”? </p> <p> Inside the ‘src_decomposition.f90’ inside the ‘ <span class="caps"> SUBROUTINE </span> decompose_cm’, there are local variables ‘nzbounds_e, nzbounds_w, nzbounds_s, nzbounds_n’. Normally these have the value 2, and cause 2 additional pixels on each side of the respective coarse grid (which are needed later on during interpolation). </p> <p> But if ‘east/…/north_add_in’ are used these ‘nzbounds_?’ are set to zero and cause (during the end of the 3rd Section of the subroutine) a too small coarse “sub”-grid. </p> <p> I changed it so that the 2 additional pixels are still taken but afterwards a check is made, that caps the value to (1,ie_in_tot) and (1,je_in_tot) </p> !&lt; ZENTEK nzbounds_e = 2 nzbounds_w = 2 nzbounds_s = 2 nzbounds_n = 2 !ZENTEK&gt; ! looking for southern grid point DO j = 1, je_in_tot IF (latitudes_in(j) &gt;= zminlat) THEN jzstart_in = j-nzbounds_s-1 EXIT ENDIF ENDDO ! looking for northern grid point DO j = 1, je_in_tot IF (latitudes_in(j) &gt; zmaxlat) THEN jzend_in = j+nzbounds_n EXIT ENDIF ENDDO ! looking for western grid point DO i = 1, ie_in_tot IF (longitudes_in(i) &gt;= zminlon) THEN izstart_in = i-nzbounds_e-1 EXIT ENDIF ENDDO ! looking for eastern grid point DO i = 1, ie_in_tot IF (longitudes_in(i) &gt; zmaxlon) THEN izend_in = i+nzbounds_w EXIT ENDIF ENDDO !&lt; ZENTEK IF (jzstart_in &lt; 1) THEN jzstart_in = 1 ENDIF IF (jzend_in &gt; ie_in_tot ) THEN jzend_in = ie_in_tot ENDIF IF (izstart_in &lt; 1) THEN izstart_in = 1 ENDIF IF (izend_in &gt; ie_in_tot) THEN izend_in = ie_in_tot ENDIF !ZENTEK&gt; <p> All in all I think these 4+3*4 lines should solve the problem. </p> <p> Leaving only the small error directly over the south pole, which can be ‘workarounded’(*) by the 3 lines of code from my last post. </p>

ok; (I think )I found the error. I would be grateful if somebody familiar with the source code could check it. Or should I contact the “Current Code Owner”?

Inside the ‘src_decomposition.f90’ inside the ‘ SUBROUTINE decompose_cm’, there are local variables ‘nzbounds_e, nzbounds_w, nzbounds_s, nzbounds_n’. Normally these have the value 2, and cause 2 additional pixels on each side of the respective coarse grid (which are needed later on during interpolation).

But if ‘east/…/north_add_in’ are used these ‘nzbounds_?’ are set to zero and cause (during the end of the 3rd Section of the subroutine) a too small coarse “sub”-grid.

I changed it so that the 2 additional pixels are still taken but afterwards a check is made, that caps the value to (1,ie_in_tot) and (1,je_in_tot)

!< ZENTEK nzbounds_e = 2 nzbounds_w = 2 nzbounds_s = 2 nzbounds_n = 2 !ZENTEK> ! looking for southern grid point DO j = 1, je_in_tot IF (latitudes_in(j) >= zminlat) THEN jzstart_in = j-nzbounds_s-1 EXIT ENDIF ENDDO ! looking for northern grid point DO j = 1, je_in_tot IF (latitudes_in(j) > zmaxlat) THEN jzend_in = j+nzbounds_n EXIT ENDIF ENDDO ! looking for western grid point DO i = 1, ie_in_tot IF (longitudes_in(i) >= zminlon) THEN izstart_in = i-nzbounds_e-1 EXIT ENDIF ENDDO ! looking for eastern grid point DO i = 1, ie_in_tot IF (longitudes_in(i) > zmaxlon) THEN izend_in = i+nzbounds_w EXIT ENDIF ENDDO !< ZENTEK IF (jzstart_in < 1) THEN jzstart_in = 1 ENDIF IF (jzend_in > ie_in_tot ) THEN jzend_in = ie_in_tot ENDIF IF (izstart_in < 1) THEN izstart_in = 1 ENDIF IF (izend_in > ie_in_tot) THEN izend_in = ie_in_tot ENDIF !ZENTEK>

All in all I think these 4+3*4 lines should solve the problem.

Leaving only the small error directly over the south pole, which can be ‘workarounded’(*) by the 3 lines of code from my last post.

<p> Thanks a lot for diving into the code B) <br/> I will check this within the next couple of days. </p>

  @burkhardtrockel in #67a0aff

<p> Thanks a lot for diving into the code B) <br/> I will check this within the next couple of days. </p>

Thanks a lot for diving into the code B)
I will check this within the next couple of days.

<p> I just tested your Pacific domain and was able to reproduce the stripes you got. The changes you made make sense. However, it should read <code> je_in_tot </code> in your listing for <code> jzend_in </code> . <br/> Since the changes do not have any effect on all of the standard tests for int2lm, I will include them in the next int2lm_clm version. <br/> Thanks for reporting and even more solving this issue! </p>

  @burkhardtrockel in #9bd1c8f

<p> I just tested your Pacific domain and was able to reproduce the stripes you got. The changes you made make sense. However, it should read <code> je_in_tot </code> in your listing for <code> jzend_in </code> . <br/> Since the changes do not have any effect on all of the standard tests for int2lm, I will include them in the next int2lm_clm version. <br/> Thanks for reporting and even more solving this issue! </p>

I just tested your Pacific domain and was able to reproduce the stripes you got. The changes you made make sense. However, it should read je_in_tot in your listing for jzend_in .
Since the changes do not have any effect on all of the standard tests for int2lm, I will include them in the next int2lm_clm version.
Thanks for reporting and even more solving this issue!

<p> Thank you for the checking and finding the <code> je_in_tot </code> -error. </p>

  @rolfzentek in #67ce79a

<p> Thank you for the checking and finding the <code> je_in_tot </code> -error. </p>

Thank you for the checking and finding the je_in_tot -error.