#2 RK3 time-stepping currently uses forcing read at n+1/2 during a time-step

src/OCE/ISF/isfhdiv.F90

@@ -71,6 +71,7 @@ CONTAINS
       !
    END SUBROUTINE isf_hdiv
 
+
    SUBROUTINE isf_hdiv_mlt(ktop, kbot, phtbl, pfrac, pfwf, pfwf_b, phdiv)
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE sbc_isf_div  ***
@@ -97,7 +98,11 @@ CONTAINS
       !
       ! compute integrated divergence correction
       DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
+#if defined key_RK3
+         zhdiv(ji,jj) = pfwf(ji,jj) * r1_rho0 / phtbl(ji,jj)
+#else
          zhdiv(ji,jj) = 0.5_wp * ( pfwf(ji,jj) + pfwf_b(ji,jj) ) * r1_rho0 / phtbl(ji,jj)
+#endif
       END_2D
       !
       ! update divergence at each level affected by ice shelf top boundary layer

src/OCE/SBC/sbcrnf.F90

@@ -212,7 +212,11 @@ CONTAINS
          IF( ln_linssh ) THEN    !* constant volume case : just apply the runoff input flow
             DO_2D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
                DO jk = 1, nk_rnf(ji,jj)
-                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact / h_rnf(ji,jj)
+#if defined key_RK3
+                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - rnf(ji,jj) * r1_rho0 / h_rnf(ji,jj)                    ! RK3: rnf forcing at n+1/2
+#else
+                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact / h_rnf(ji,jj)   ! MLF: rnf forcing at Kmm (n)
+#endif
                END DO
             END_2D
          ELSE                    !* variable volume case
@@ -224,7 +228,11 @@ CONTAINS
             END_2D
             DO_2D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls )         ! apply the runoff input flow
                DO jk = 1, nk_rnf(ji,jj)
-                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact / h_rnf(ji,jj)
+#if defined key_RK3
+                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - rnf(ji,jj) * r1_rho0 / h_rnf(ji,jj)                    ! RK3: rnf forcing at n+1/2
+#else
+                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact / h_rnf(ji,jj)   ! MLF: rnf forcing at Kmm (n)
+#endif
                END DO
             END_2D
          ENDIF
@@ -233,7 +241,11 @@ CONTAINS
             h_rnf (ji,jj)   = e3t(ji,jj,1,Kmm)        ! update h_rnf to be depth of top box
          END_2D
          DO_2D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
-            phdivn(ji,jj,1) = phdivn(ji,jj,1) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact / e3t(ji,jj,1,Kmm)
+#if defined key_RK3
+            phdivn(ji,jj,1) = phdivn(ji,jj,1) - rnf(ji,jj) * r1_rho0 / e3t(ji,jj,1,Kmm)                    ! RK3: rnf forcing at n+1/2
+#else
+            phdivn(ji,jj,1) = phdivn(ji,jj,1) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact / e3t(ji,jj,1,Kmm)   ! MLF: rnf forcing at Kmm (n)
+#endif
          END_2D
       ENDIF
       !

src/OCE/SBC/sbcssm.F90

@@ -79,7 +79,11 @@ CONTAINS
          ENDIF
          sss_m(:,:) = zts(:,:,jp_sal)
          !                          ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
-         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = ssh(:,:,Kmm) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
+#if defined key_RK3
+         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = ssh(:,:,Kmm) - ssh_ib(:,:)                            ! RK3: forcing at n+1/2
+#else
+         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = ssh(:,:,Kmm) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )   ! MLF: forcing at n (Kmm)
+#endif
          ELSE                    ;   ssh_m(:,:) = ssh(:,:,Kmm)
          ENDIF
          !
@@ -102,7 +106,11 @@ CONTAINS
             ENDIF
             sss_m(:,:) = zcoef * zts(:,:,jp_sal)
             !                          ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
+#if defined key_RK3
+            IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = zcoef * ( ssh(:,:,Kmm) - ssh_ib(:,:) )
+#else
             IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = zcoef * ( ssh(:,:,Kmm) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) )
+#endif
             ELSE                    ;   ssh_m(:,:) = zcoef *   ssh(:,:,Kmm)
             ENDIF
             !
@@ -130,7 +138,11 @@ CONTAINS
          ENDIF
          sss_m(:,:) = sss_m(:,:) + zts(:,:,jp_sal)
          !                          ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
+#if defined key_RK3
+         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = ssh_m(:,:) + ssh(:,:,Kmm) - ssh_ib(:,:)
+#else
          IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = ssh_m(:,:) + ssh(:,:,Kmm) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
+#endif
          ELSE                    ;   ssh_m(:,:) = ssh_m(:,:) + ssh(:,:,Kmm)
          ENDIF
          !