Investigate halo1 - halo2 differences with ICB

[acc]

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:

icbutl.F90

@@ -279 +365 @@
          ! note that here there is no +0.5 added
          ! since we're looking for four T points containing quadrant we're in of
          ! current T cell
-         kii = MAX(0, INT( pi        ))
-         kij = MAX(0, INT( pj        ))    ! T-point
-         zwi = pi - REAL(kii,wp)
-         zwj = pj - REAL(kij,wp)
+         kii = MAX(0, INT( pi-(nn_hls-1)        ))
+         kij = MAX(0, INT( pj-(nn_hls-1)        ))    ! T-point
+         zwi = pi-(nn_hls-1) - REAL(kii,wp)
+         zwj = pj-(nn_hls-1) - REAL(kij,wp)
       CASE ( 'U' )
-         kii = MAX(0, INT( pi-0.5_wp ))
-         kij = MAX(0, INT( pj        ))    ! U-point
-         zwi = pi - 0.5_wp - REAL(kii,wp)
-         zwj = pj - REAL(kij,wp)
+         kii = MAX(0, INT( pi-(nn_hls-1)-0.5_wp ))
+         kij = MAX(0, INT( pj-(nn_hls-1)        ))    ! U-point
+         zwi = pi-(nn_hls-1) - 0.5_wp - REAL(kii,wp)
+         zwj = pj-(nn_hls-1) - REAL(kij,wp)
       CASE ( 'V' )
-         kii = MAX(0, INT( pi        ))
-         kij = MAX(0, INT( pj-0.5_wp ))    ! V-point
-         zwi = pi - REAL(kii,wp)
-         zwj = pj - 0.5_wp - REAL(kij,wp)
+         kii = MAX(0, INT( pi-(nn_hls-1)        ))
+         kij = MAX(0, INT( pj-(nn_hls-1)-0.5_wp ))    ! V-point
+         zwi = pi-(nn_hls-1) - REAL(kii,wp)
+         zwj = pj-(nn_hls-1) - 0.5_wp - REAL(kij,wp)
       CASE ( 'F' )
-         kii = MAX(0, INT( pi-0.5_wp ))
-         kij = MAX(0, INT( pj-0.5_wp ))    ! F-point
-         zwi = pi - 0.5_wp - REAL(kii,wp)
-         zwj = pj - 0.5_wp - REAL(kij,wp)
+         kii = MAX(0, INT( pi-(nn_hls-1)-0.5_wp ))
+         kij = MAX(0, INT( pj-(nn_hls-1)-0.5_wp ))    ! F-point
+         zwi = pi-(nn_hls-1) - 0.5_wp - REAL(kii,wp)
+         zwj = pj-(nn_hls-1) - 0.5_wp - REAL(kij,wp)
       END SELECT
+      kii = kii + (nn_hls-1)
+      kij = kij + (nn_hls-1)
       !
       ! compute weight
       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.:

icbdyn.F90

@@ -96 +96 @@
          ! position using di/dt & djdt   !   V2  in m/s
          zxi2 = zxi1 + zdt_2 * zu1          ;   zuvel2 = zuvel1 + zdt_2 * zax1
          zyj2 = zyj1 + zdt_2 * zv1          ;   zvvel2 = zvvel1 + zdt_2 * zay1
+         zuvel2 = int(zuvel2 * 1.e10_wp,KIND=8)*1.e-10_wp
+         zvvel2 = int(zvvel2 * 1.e10_wp,KIND=8)*1.e-10_wp
+         zxi2 = int(zxi2 * 1.e10_wp,KIND=8)*1.e-10_wp
+         zyj2 = int(zyj2 * 1.e10_wp,KIND=8)*1.e-10_wp
          !
          CALL icb_ground( berg, zxi2, zxi1, zu1,   &
             &                   zyj2, zyj1, zv1, ll_bounced )
@@ -112 +116 @@
          !                                         !**  X3 = X1+dt/2*V2  ;   V3 = V1+dt/2*A2; A3=A(X3)
          zxi3  = zxi1  + zdt_2 * zu2   ;   zuvel3 = zuvel1 + zdt_2 * zax2
          zyj3  = zyj1  + zdt_2 * zv2   ;   zvvel3 = zvvel1 + zdt_2 * zay2
+         zuvel3 = int(zuvel3 * 1.e10_wp,KIND=8)*1.e-10_wp
+         zvvel3 = int(zvvel3 * 1.e10_wp,KIND=8)*1.e-10_wp
+         zxi3 = int(zxi3 * 1.e10_wp,KIND=8)*1.e-10_wp
+         zyj3 = int(zyj3 * 1.e10_wp,KIND=8)*1.e-10_wp
          !
          CALL icb_ground( berg, zxi3, zxi1, zu2,   &
             &                   zyj3, zyj1, zv2, ll_bounced )
@@ -128 +136 @@
          !                                         !**   X4 = X1+dt*V3   ;   V4 = V1+dt*A3
          zxi4 = zxi1 + zdt * zu3   ;   zuvel4 = zuvel1 + zdt * zax3
          zyj4 = zyj1 + zdt * zv3   ;   zvvel4 = zvvel1 + zdt * zay3
+         zuvel4 = int(zuvel4 * 1.e10_wp,KIND=8)*1.e-10_wp
+         zvvel4 = int(zvvel4 * 1.e10_wp,KIND=8)*1.e-10_wp
+         zxi4 = int(zxi4 * 1.e10_wp,KIND=8)*1.e-10_wp
+         zyj4 = int(zyj4 * 1.e10_wp,KIND=8)*1.e-10_wp
 
          CALL icb_ground( berg, zxi4, zxi1, zu3,   &
             &                   zyj4, zyj1, zv3, ll_bounced )
@@ -151 +163 @@
          CALL icb_ground( berg, zxi_n, zxi1, zuvel_n,   &
             &                   zyj_n, zyj1, zvvel_n, ll_bounced )
 
-         pt%uvel = zuvel_n                        !** save in berg structure
-         pt%vvel = zvvel_n
-         pt%xi   = zxi_n
-         pt%yj   = zyj_n
+         !pt%uvel = zuvel_n                        !** save in berg structure
+         !pt%vvel = zvvel_n
+         pt%uvel = int(zuvel_n * 1.e10_wp,KIND=8)*1.e-10_wp                        !** save in berg structure
+         pt%vvel = int(zvvel_n * 1.e10_wp,KIND=8)*1.e-10_wp
+         pt%xi   = int( zxi_n * 1.e10_wp,KIND=8)*1.e-10_wp
+         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.

[acc]

In 15040:

Error: Failed to load processor CommitTicketReference

No macro or processor named 'CommitTicketReference' found

[acc]

In 15043:

Error: Failed to load processor CommitTicketReference

No macro or processor named 'CommitTicketReference' found

[acc]

In 15088:

Error: Failed to load processor CommitTicketReference

No macro or processor named 'CommitTicketReference' found