| Description |
Context
SETTE tests with ORCA2_ICE_PISCES give different results with nn_hls=1 and nn_hls=2 if icebergs are active. This ticket is to record investigations and potential solutions.
Analysis
This is a tricky one to pin down because differences in run.stat outputs don't show for many timesteps (84 on my test system), the initial differences are very small and much of the calculation in icb routines is done through function calls that are difficult to instrument. First observations are:
- Running with icebergs on but ln_passive_mode =.true. will give identical results. This isn't surprising since the icebergs do not affect the ocean in this mode. It does, however, provide an option for running with icebergs on for the purposes of assessing computational performance with nn_hls=2 vs nn_hls=1.
- With ddefault settings, the earliest detectable diffferences can be found after 6 or 7 timesteps in the values returned from icb_accel (icbdyn.F90). These are caused by 14th or 15th decimal place differences in the interpolation weights returned by icb_utl_pos which is called from icb_utl_interp which is called by icb_accel. These differences occur despite the iceberg positions appearing to be 'identical' in the two cases. iceberg positions are given as a continuous local gridcell value so they are offset by one between the two cases but appear otherwise identical to 15 decimal places. Nevertheless the interpolation weights in the two cases have small differences which eventually lead to significant differences. The onset of differences can be delayed by ensuring the calculation of weights is identical in both cases by removing the additional offset in the nn_hls=2 case. For example:
-
svn diff icbutl.F90
|
|
|
|
| 279 | 365 | ! note that here there is no +0.5 added |
| 280 | 366 | ! since we're looking for four T points containing quadrant we're in of |
| 281 | 367 | ! current T cell |
| 282 | | kii = MAX(0, INT( pi |
| 283 | | kij = MAX(0, INT( pj |
| 284 | | zwi = pi |
| 285 | | zwj = pj |
| | 368 | kii = MAX(0, INT( pi-(nn_hls-1) )) |
| | 369 | kij = MAX(0, INT( pj-(nn_hls-1) )) ! T-point |
| | 370 | zwi = pi-(nn_hls-1) - REAL(kii,wp) |
| | 371 | zwj = pj-(nn_hls-1) - REAL(kij,wp) |
| 286 | 372 | CASE ( 'U' ) |
| 287 | | kii = MAX(0, INT( pi- |
| 288 | | kij = MAX(0, INT( pj |
| 289 | | zwi = pi |
| 290 | | zwj = pj |
| | 373 | kii = MAX(0, INT( pi-(nn_hls-1)-0.5_wp )) |
| | 374 | kij = MAX(0, INT( pj-(nn_hls-1) )) ! U-point |
| | 375 | zwi = pi-(nn_hls-1) - 0.5_wp - REAL(kii,wp) |
| | 376 | zwj = pj-(nn_hls-1) - REAL(kij,wp) |
| 291 | 377 | CASE ( 'V' ) |
| 292 | | kii = MAX(0, INT( pi |
| 293 | | kij = MAX(0, INT( pj- |
| 294 | | zwi = pi |
| 295 | | zwj = pj |
| | 378 | kii = MAX(0, INT( pi-(nn_hls-1) )) |
| | 379 | kij = MAX(0, INT( pj-(nn_hls-1)-0.5_wp )) ! V-point |
| | 380 | zwi = pi-(nn_hls-1) - REAL(kii,wp) |
| | 381 | zwj = pj-(nn_hls-1) - 0.5_wp - REAL(kij,wp) |
| 296 | 382 | CASE ( 'F' ) |
| 297 | | kii = MAX(0, INT( pi- |
| 298 | | kij = MAX(0, INT( pj- |
| 299 | | zwi = pi |
| 300 | | zwj = pj |
| | 383 | kii = MAX(0, INT( pi-(nn_hls-1)-0.5_wp )) |
| | 384 | kij = MAX(0, INT( pj-(nn_hls-1)-0.5_wp )) ! F-point |
| | 385 | zwi = pi-(nn_hls-1) - 0.5_wp - REAL(kii,wp) |
| | 386 | zwj = pj-(nn_hls-1) - 0.5_wp - REAL(kij,wp) |
| 301 | 387 | END SELECT |
| | 388 | kii = kii + (nn_hls-1) |
| | 389 | kij = kij + (nn_hls-1) |
| 302 | 390 | ! |
| 303 | 391 | ! compute weight |
| 304 | 392 | pw(1) = (1._wp-zwi) * (1._wp-zwj) |
There is further complication in that icb_utl_interp also estimates the local ssh gradient using:
CALL icb_utl_pos( pi+0.1_wp, pj, 'T', iiTp, ijTp, zwTp, zmskTp )
CALL icb_utl_pos( pi-0.1_wp, pj, 'T', iiTm, ijTm, zwTm, zmskTm )
!
IF ( .NOT. PRESENT(pe1) ) pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj )
pssh_i = ( icb_utl_bilin_h( ssh_e, iiTp, ijTp, zwTp*zmskTp, .false. ) - &
& icb_utl_bilin_h( ssh_e, iiTm, ijTm, zwTm*zmskTm, .false. ) ) / ( 0.2_wp * pe1 )
To make these calculations identical, the (nn_hls-1) factor has to be subtracted before the offset is applied. An alternative version of icb_utl_pos can be created for this purpose.
These changes further delay the onset of differences but do not eliminate them.
- A test has been carried out where the least signifcant decimal places of the iceberg positions and speeds are set to zero. I.e.:
-
svn diff icbdyn.F90
|
|
|
|
| 96 | 96 | ! position using di/dt & djdt ! V2 in m/s |
| 97 | 97 | zxi2 = zxi1 + zdt_2 * zu1 ; zuvel2 = zuvel1 + zdt_2 * zax1 |
| 98 | 98 | zyj2 = zyj1 + zdt_2 * zv1 ; zvvel2 = zvvel1 + zdt_2 * zay1 |
| | 99 | zuvel2 = int(zuvel2 * 1.e10_wp,KIND=8)*1.e-10_wp |
| | 100 | zvvel2 = int(zvvel2 * 1.e10_wp,KIND=8)*1.e-10_wp |
| | 101 | zxi2 = int(zxi2 * 1.e10_wp,KIND=8)*1.e-10_wp |
| | 102 | zyj2 = int(zyj2 * 1.e10_wp,KIND=8)*1.e-10_wp |
| 99 | 103 | ! |
| 100 | 104 | CALL icb_ground( berg, zxi2, zxi1, zu1, & |
| 101 | 105 | & zyj2, zyj1, zv1, ll_bounced ) |
| … |
… |
|
| 112 | 116 | ! !** X3 = X1+dt/2*V2 ; V3 = V1+dt/2*A2; A3=A(X3) |
| 113 | 117 | zxi3 = zxi1 + zdt_2 * zu2 ; zuvel3 = zuvel1 + zdt_2 * zax2 |
| 114 | 118 | zyj3 = zyj1 + zdt_2 * zv2 ; zvvel3 = zvvel1 + zdt_2 * zay2 |
| | 119 | zuvel3 = int(zuvel3 * 1.e10_wp,KIND=8)*1.e-10_wp |
| | 120 | zvvel3 = int(zvvel3 * 1.e10_wp,KIND=8)*1.e-10_wp |
| | 121 | zxi3 = int(zxi3 * 1.e10_wp,KIND=8)*1.e-10_wp |
| | 122 | zyj3 = int(zyj3 * 1.e10_wp,KIND=8)*1.e-10_wp |
| 115 | 123 | ! |
| 116 | 124 | CALL icb_ground( berg, zxi3, zxi1, zu2, & |
| 117 | 125 | & zyj3, zyj1, zv2, ll_bounced ) |
| … |
… |
|
| 128 | 136 | ! !** X4 = X1+dt*V3 ; V4 = V1+dt*A3 |
| 129 | 137 | zxi4 = zxi1 + zdt * zu3 ; zuvel4 = zuvel1 + zdt * zax3 |
| 130 | 138 | zyj4 = zyj1 + zdt * zv3 ; zvvel4 = zvvel1 + zdt * zay3 |
| | 139 | zuvel4 = int(zuvel4 * 1.e10_wp,KIND=8)*1.e-10_wp |
| | 140 | zvvel4 = int(zvvel4 * 1.e10_wp,KIND=8)*1.e-10_wp |
| | 141 | zxi4 = int(zxi4 * 1.e10_wp,KIND=8)*1.e-10_wp |
| | 142 | zyj4 = int(zyj4 * 1.e10_wp,KIND=8)*1.e-10_wp |
| 131 | 143 | |
| 132 | 144 | CALL icb_ground( berg, zxi4, zxi1, zu3, & |
| 133 | 145 | & zyj4, zyj1, zv3, ll_bounced ) |
| … |
… |
|
| 151 | 163 | CALL icb_ground( berg, zxi_n, zxi1, zuvel_n, & |
| 152 | 164 | & zyj_n, zyj1, zvvel_n, ll_bounced ) |
| 153 | 165 | |
| 154 | | pt%uvel = zuvel_n !** save in berg structure |
| 155 | | pt%vvel = zvvel_n |
| 156 | | pt%xi = zxi_n |
| 157 | | pt%yj = zyj_n |
| | 166 | !pt%uvel = zuvel_n !** save in berg structure |
| | 167 | !pt%vvel = zvvel_n |
| | 168 | pt%uvel = int(zuvel_n * 1.e10_wp,KIND=8)*1.e-10_wp !** save in berg structure |
| | 169 | pt%vvel = int(zvvel_n * 1.e10_wp,KIND=8)*1.e-10_wp |
| | 170 | pt%xi = int( zxi_n * 1.e10_wp,KIND=8)*1.e-10_wp |
| | 171 | pt%yj = int( zyj_n * 1.e10_wp,KIND=8)*1.e-10_wp |
With these changes differences in run.stat values do not occur until timestep 255.
Recommendation
It seems the differences are due to numerical differences arising from the offset in local gridcell number. A solution is to remove the factor of (nn_hls-1) from the initial iceberg positions and to add this factor back on whenever continuous positions are converted back to gridcell indices. Doing this ensures that all positional calculations are done with exactly the same numbers and results in identical solutions with either nn_hls=1 or nn_hls=2. Changes have been made and tested on a ticket branch pending review and independent testing.
|
