Index: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynspg_ts.F90	(revision 13650)
+++ NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynspg_ts.F90	(revision 14017)
@@ -1,5 +1,5 @@
 MODULE dynspg_ts
 
-   !! Includes ROMS wd scheme with diagnostic outputs ; puu(:,:,:,Kmm) and puu(:,:,:,Krhs) updates are commented out ! 
+   !! Includes ROMS wd scheme with diagnostic outputs ; puu(:,:,:,Kmm) and puu(:,:,:,Krhs) updates are commented out !
 
    !!======================================================================
@@ -23,5 +23,5 @@
 
    !!----------------------------------------------------------------------
-   !!   dyn_spg_ts     : compute surface pressure gradient trend using a time-splitting scheme 
+   !!   dyn_spg_ts     : compute surface pressure gradient trend using a time-splitting scheme
    !!   dyn_spg_ts_init: initialisation of the time-splitting scheme
    !!   ts_wgt         : set time-splitting weights for temporal averaging (or not)
@@ -50,5 +50,5 @@
    USE agrif_oce
 #endif
-#if defined key_asminc   
+#if defined key_asminc
    USE asminc          ! Assimilation increment
 #endif
@@ -66,5 +66,5 @@
    PRIVATE
 
-   PUBLIC dyn_spg_ts        ! called by dyn_spg 
+   PUBLIC dyn_spg_ts        ! called by dyn_spg
    PUBLIC dyn_spg_ts_init   !    -    - dyn_spg_init
 
@@ -122,10 +122,10 @@
       !! ** Purpose : - Compute the now trend due to the explicit time stepping
       !!              of the quasi-linear barotropic system, and add it to the
-      !!              general momentum trend. 
+      !!              general momentum trend.
       !!
       !! ** Method  : - split-explicit schem (time splitting) :
       !!      Barotropic variables are advanced from internal time steps
       !!      "n"   to "n+1" if ln_bt_fw=T
-      !!      or from 
+      !!      or from
       !!      "n-1" to "n+1" if ln_bt_fw=F
       !!      thanks to a generalized forward-backward time stepping (see ref. below).
@@ -136,9 +136,9 @@
       !!      -Compute barotropic advective fluxes at step "n"     : un_adv, vn_adv
       !!      These are used to advect tracers and are compliant with discrete
-      !!      continuity equation taken at the baroclinic time steps. This 
+      !!      continuity equation taken at the baroclinic time steps. This
       !!      ensures tracers conservation.
       !!      - (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) momentum trend updated with barotropic component.
       !!
-      !! References : Shchepetkin and McWilliams, Ocean Modelling, 2005. 
+      !! References : Shchepetkin and McWilliams, Ocean Modelling, 2005.
       !!---------------------------------------------------------------------
       INTEGER                             , INTENT( in )  ::  kt                  ! ocean time-step index
@@ -148,5 +148,5 @@
       !
       INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices
-      LOGICAL  ::   ll_fw_start           ! =T : forward integration 
+      LOGICAL  ::   ll_fw_start           ! =T : forward integration
       LOGICAL  ::   ll_init               ! =T : special startup of 2d equations
       INTEGER  ::   noffset               ! local integers  : time offset for bdy update
@@ -164,5 +164,5 @@
       REAL(wp), DIMENSION(jpi, jpj, jpk) :: ze3u, ze3v
       !
-      REAL(wp) ::   zwdramp                     ! local scalar - only used if ln_wd_dl = .True. 
+      REAL(wp) ::   zwdramp                     ! local scalar - only used if ln_wd_dl = .True.
 
       INTEGER  :: iwdg, jwdg, kwdg   ! short-hand values for the indices of the output point
@@ -179,14 +179,14 @@
       IF( ln_wd_dl ) ALLOCATE( ztwdmask(jpi,jpj), zuwdmask(jpi,jpj), zvwdmask(jpi,jpj), zuwdav2(jpi,jpj), zvwdav2(jpi,jpj))
       !
-      zwdramp = r_rn_wdmin1               ! simplest ramp 
+      zwdramp = r_rn_wdmin1               ! simplest ramp
 !     zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp
-      !                                         ! inverse of baroclinic time step 
-      r1_Dt_b = 1._wp / rDt 
-      !
-      ll_init     = ln_bt_av                    ! if no time averaging, then no specific restart 
+      !                                         ! inverse of baroclinic time step
+      r1_Dt_b = 1._wp / rDt
+      !
+      ll_init     = ln_bt_av                    ! if no time averaging, then no specific restart
       ll_fw_start = .FALSE.
       !                                         ! time offset in steps for bdy data update
       IF( .NOT.ln_bt_fw ) THEN   ;   noffset = - nn_e
-      ELSE                       ;   noffset =   0 
+      ELSE                       ;   noffset =   0
       ENDIF
       !
@@ -215,5 +215,5 @@
             ! and the averaging weights. We don't have an easy way of telling whether we did
             ! an Euler timestep on the first timestep (because l_1st_euler is reset to .false.
-            ! at the end of the first timestep) so just do this in all cases. 
+            ! at the end of the first timestep) so just do this in all cases.
             ll_fw_start = .FALSE.
             CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
@@ -225,5 +225,5 @@
       !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
       ! -----------------------------------------------------------------------------
-      !      
+      !
       !
       !                                   !=  zu_frc =  1/H e3*d/dt(Ua)  =!  (Vertical mean of Ua, the 3D trends)
@@ -243,8 +243,8 @@
          vv(:,:,jk,Krhs) = ( vv(:,:,jk,Krhs) - zv_frc(:,:) ) * vmask(:,:,jk)
       END DO
-      
+
 !!gm  Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum....
 !!gm  Is it correct to do so ?   I think so...
-      
+
       !                                   !=  remove 2D Coriolis and pressure gradient trends  =!
       !                                   !  -------------------------------------------------  !
@@ -254,5 +254,5 @@
       !
       !                                         !* 2D Coriolis trends
-      zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:)        ! now fluxes 
+      zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:)        ! now fluxes
       zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:)        ! NB: FULL domain : put a value in last row and column
       !
@@ -273,5 +273,5 @@
             DO_2D( 0, 0, 0, 0 )
                zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  pssh(ji+1,jj  ,Kmm) - pssh(ji  ,jj  ,Kmm)  ) * r1_e1u(ji,jj)
-               zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  pssh(ji  ,jj+1,Kmm) - pssh(ji  ,jj  ,Kmm)  ) * r1_e2v(ji,jj) 
+               zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  pssh(ji  ,jj+1,Kmm) - pssh(ji  ,jj  ,Kmm)  ) * r1_e2v(ji,jj)
             END_2D
          ENDIF
@@ -319,5 +319,5 @@
             zv_frc(ji,jj) =  zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kmm)
          END_2D
-      ENDIF  
+      ENDIF
       !
       !              !----------------!
@@ -369,11 +369,11 @@
       !
       ! -----------------------------------------------------------------------
-      !  Phase 2 : Integration of the barotropic equations 
+      !  Phase 2 : Integration of the barotropic equations
       ! -----------------------------------------------------------------------
       !
       !                                             ! ==================== !
       !                                             !    Initialisations   !
-      !                                             ! ==================== !  
-      ! Initialize barotropic variables:      
+      !                                             ! ==================== !
+      ! Initialize barotropic variables:
       IF( ll_init )THEN
          sshbb_e(:,:) = 0._wp
@@ -391,26 +391,26 @@
       ENDIF
       !
-      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                    
-         sshn_e(:,:) =    pssh(:,:,Kmm)            
-         un_e  (:,:) =    puu_b(:,:,Kmm)            
+      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields
+         sshn_e(:,:) =    pssh(:,:,Kmm)
+         un_e  (:,:) =    puu_b(:,:,Kmm)
          vn_e  (:,:) =    pvv_b(:,:,Kmm)
          !
-         hu_e  (:,:) =    hu(:,:,Kmm)       
-         hv_e  (:,:) =    hv(:,:,Kmm) 
-         hur_e (:,:) = r1_hu(:,:,Kmm)    
+         hu_e  (:,:) =    hu(:,:,Kmm)
+         hv_e  (:,:) =    hv(:,:,Kmm)
+         hur_e (:,:) = r1_hu(:,:,Kmm)
          hvr_e (:,:) = r1_hv(:,:,Kmm)
       ELSE                                ! CENTRED integration: start from BEFORE fields
          sshn_e(:,:) =    pssh(:,:,Kbb)
-         un_e  (:,:) =    puu_b(:,:,Kbb)         
+         un_e  (:,:) =    puu_b(:,:,Kbb)
          vn_e  (:,:) =    pvv_b(:,:,Kbb)
          !
-         hu_e  (:,:) =    hu(:,:,Kbb)       
-         hv_e  (:,:) =    hv(:,:,Kbb) 
-         hur_e (:,:) = r1_hu(:,:,Kbb)    
+         hu_e  (:,:) =    hu(:,:,Kbb)
+         hv_e  (:,:) =    hv(:,:,Kbb)
+         hur_e (:,:) = r1_hu(:,:,Kbb)
          hvr_e (:,:) = r1_hv(:,:,Kbb)
       ENDIF
       !
       ! Initialize sums:
-      puu_b  (:,:,Kaa) = 0._wp       ! After barotropic velocities (or transport if flux form)          
+      puu_b  (:,:,Kaa) = 0._wp       ! After barotropic velocities (or transport if flux form)
       pvv_b  (:,:,Kaa) = 0._wp
       pssh  (:,:,Kaa) = 0._wp       ! Sum for after averaged sea level
@@ -419,9 +419,9 @@
       !
       IF( ln_wd_dl ) THEN
-         zuwdmask(:,:) = 0._wp  ! set to zero for definiteness (not sure this is necessary) 
-         zvwdmask(:,:) = 0._wp  ! 
-         zuwdav2 (:,:) = 0._wp 
-         zvwdav2 (:,:) = 0._wp   
-      END IF 
+         zuwdmask(:,:) = 0._wp  ! set to zero for definiteness (not sure this is necessary)
+         zvwdmask(:,:) = 0._wp  !
+         zuwdav2 (:,:) = 0._wp
+         zvwdav2 (:,:) = 0._wp
+      END IF
 
       !                                             ! ==================== !
@@ -443,9 +443,9 @@
          !
          !                       !* Set extrapolation coefficients for predictor step:
-         IF ((jn<3).AND.ll_init) THEN      ! Forward           
-           za1 = 1._wp                                          
-           za2 = 0._wp                        
-           za3 = 0._wp                        
-         ELSE                              ! AB3-AM4 Coefficients: bet=0.281105 
+         IF ((jn<3).AND.ll_init) THEN      ! Forward
+           za1 = 1._wp
+           za2 = 0._wp
+           za3 = 0._wp
+         ELSE                              ! AB3-AM4 Coefficients: bet=0.281105
            za1 =  1.781105_wp              ! za1 =   3/2 +   bet
            za2 = -1.06221_wp               ! za2 = -(1/2 + 2*bet)
@@ -467,5 +467,5 @@
             !--------------------------------------------------------------------------------!
             zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
-            
+
             ! set wetting & drying mask at tracer points for this barotropic mid-step
             IF( ln_wd_dl )   CALL wad_tmsk( zsshp2_e, ztwdmask )
@@ -493,5 +493,5 @@
          !                             ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d
          IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e )
-         !      
+         !
          !                             ! resulting flux at mid-step (not over the full domain)
          zhU(1:jpim1,1:jpj  ) = e2u(1:jpim1,1:jpj  ) * ua_e(1:jpim1,1:jpj  ) * zhup2_e(1:jpim1,1:jpj  )   ! not jpi-column
@@ -504,9 +504,9 @@
          IF( ln_wd_il )   CALL wad_lmt_bt(zhU, zhV, sshn_e, zssh_frc, rDt_e)    !!gm wad_lmt_bt use of lbc_lnk on zhU, zhV
 
-         IF( ln_wd_dl ) THEN           ! un_e and vn_e are set to zero at faces where 
+         IF( ln_wd_dl ) THEN           ! un_e and vn_e are set to zero at faces where
             !                          ! the direction of the flow is from dry cells
             CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask )   ! not jpi colomn for U, not jpj row for V
             !
-         ENDIF    
+         ENDIF
          !
          !
@@ -532,5 +532,5 @@
          un_adv(:,:) = un_adv(:,:) + za2 * zhU(:,:) * r1_e2u(:,:)
          vn_adv(:,:) = vn_adv(:,:) + za2 * zhV(:,:) * r1_e1v(:,:)
-         ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True) 
+         ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True)
          IF ( ln_wd_dl_bc ) THEN
             zuwdav2(1:jpim1,1:jpj  ) = zuwdav2(1:jpim1,1:jpj  ) + za2 * zuwdmask(1:jpim1,1:jpj  )   ! not jpi-column
@@ -538,7 +538,7 @@
          END IF
          !
-         !  
+         !
          ! Sea Surface Height at u-,v-points (vvl case only)
-         IF( .NOT.ln_linssh ) THEN                                
+         IF( .NOT.ln_linssh ) THEN
             DO_2D( 0, 0, 0, 0 )
                zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj)    &
@@ -549,6 +549,6 @@
                   &              +   e1e2t(ji  ,jj+1)  * ssha_e(ji  ,jj+1) )
             END_2D
-         ENDIF   
-         !         
+         ENDIF
+         !
          ! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2
          !--            m+1/2           m+1              m               m-1              m-2     --!
@@ -607,8 +607,8 @@
          IF( ln_dynadv_vec .OR. ln_linssh ) THEN      !* Vector form
             DO_2D( 0, 0, 0, 0 )
-               ua_e(ji,jj) = (                                 un_e(ji,jj)   & 
+               ua_e(ji,jj) = (                                 un_e(ji,jj)   &
                          &     + rDt_e * (                   zu_spg(ji,jj)   &
                          &                                 + zu_trd(ji,jj)   &
-                         &                                 + zu_frc(ji,jj) ) & 
+                         &                                 + zu_frc(ji,jj) ) &
                          &   ) * ssumask(ji,jj)
 
@@ -632,5 +632,5 @@
                z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) )
                !
-               ua_e(ji,jj) = (               hu_e  (ji,jj) *   un_e (ji,jj)      & 
+               ua_e(ji,jj) = (               hu_e  (ji,jj) *   un_e (ji,jj)      &
                     &            + rDt_e * (  zhu_bck        * zu_spg (ji,jj)  &   !
                     &                       + zhup2_e(ji,jj) * zu_trd (ji,jj)  &   !
@@ -650,5 +650,5 @@
             END_2D
          ENDIF
-       
+
          IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
             hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1)
@@ -667,5 +667,5 @@
          !                                                 ! open boundaries
          IF( ln_bdy )   CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
-#if defined key_agrif                                                           
+#if defined key_agrif
          IF( .NOT.Agrif_Root() )  CALL agrif_dyn_ts( jn )  ! Agrif
 #endif
@@ -686,16 +686,16 @@
          !                                             !* Sum over whole bt loop
          !                                             !  ----------------------
-         za1 = wgtbtp1(jn)                                    
+         za1 = wgtbtp1(jn)
          IF( ln_dynadv_vec .OR. ln_linssh ) THEN    ! Sum velocities
-            puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) 
-            pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) 
+            puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:)
+            pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:)
          ELSE                                       ! Sum transports
-            IF ( .NOT.ln_wd_dl ) THEN  
+            IF ( .NOT.ln_wd_dl ) THEN
                puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) * hu_e (:,:)
                pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) * hv_e (:,:)
-            ELSE 
+            ELSE
                puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) * hu_e (:,:) * zuwdmask(:,:)
                pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) * hv_e (:,:) * zvwdmask(:,:)
-            END IF 
+            END IF
          ENDIF
          !                                          ! Sum sea level
@@ -715,5 +715,5 @@
                zun_save = un_adv(ji,jj)
                zvn_save = vn_adv(ji,jj)
-               !                          ! apply the previously computed correction 
+               !                          ! apply the previously computed correction
                un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - rn_atfp * un_bf(ji,jj) )
                vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - rn_atfp * vn_bf(ji,jj) )
@@ -763,5 +763,5 @@
 
 
-      ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)  
+      ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)
       DO jk = 1, jpkm1
          puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) + un_adv(:,:)*r1_hu(:,:,Kmm) - puu_b(:,:,Kmm) ) * umask(:,:,jk)
@@ -769,14 +769,14 @@
       END DO
 
-      IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN 
+      IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
          DO jk = 1, jpkm1
             puu(:,:,jk,Kmm) = ( un_adv(:,:)*r1_hu(:,:,Kmm) &
-                       & + zuwdav2(:,:)*(puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm)) ) * umask(:,:,jk) 
-            pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) & 
-                       & + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm)) ) * vmask(:,:,jk)  
+                       & + zuwdav2(:,:)*(puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm)) ) * umask(:,:,jk)
+            pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) &
+                       & + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm)) ) * vmask(:,:,jk)
          END DO
-      END IF 
-
-      
+      END IF
+
+
       CALL iom_put(  "ubar", un_adv(:,:)*r1_hu(:,:,Kmm) )    ! barotropic i-current
       CALL iom_put(  "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) )    ! barotropic i-current
@@ -796,5 +796,5 @@
          vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
       ENDIF
-#endif      
+#endif
       !                                   !* write time-spliting arrays in the restart
       IF( lrst_oce .AND.ln_bt_fw )   CALL ts_rst( kt, 'WRITE' )
@@ -817,8 +817,8 @@
       LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true.
       LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true.
-      INTEGER, INTENT(inout) :: jpit      ! cycle length    
+      INTEGER, INTENT(inout) :: jpit      ! cycle length
       REAL(wp), DIMENSION(3*nn_e), INTENT(inout) ::   zwgt1, & ! Primary weights
                                                          zwgt2    ! Secondary weights
-      
+
       INTEGER ::  jic, jn, ji                      ! temporary integers
       REAL(wp) :: za1, za2
@@ -829,5 +829,5 @@
 
       ! Set time index when averaged value is requested
-      IF (ll_fw) THEN 
+      IF (ll_fw) THEN
          jic = nn_e
       ELSE
@@ -837,6 +837,6 @@
       ! Set primary weights:
       IF (ll_av) THEN
-           ! Define simple boxcar window for primary weights 
-           ! (width = nn_e, centered around jic)     
+           ! Define simple boxcar window for primary weights
+           ! (width = nn_e, centered around jic)
          SELECT CASE ( nn_bt_flt )
               CASE( 0 )  ! No averaging
@@ -846,5 +846,5 @@
               CASE( 1 )  ! Boxcar, width = nn_e
                  DO jn = 1, 3*nn_e
-                    za1 = ABS(float(jn-jic))/float(nn_e) 
+                    za1 = ABS(float(jn-jic))/float(nn_e)
                     IF (za1 < 0.5_wp) THEN
                       zwgt1(jn) = 1._wp
@@ -855,5 +855,5 @@
               CASE( 2 )  ! Boxcar, width = 2 * nn_e
                  DO jn = 1, 3*nn_e
-                    za1 = ABS(float(jn-jic))/float(nn_e) 
+                    za1 = ABS(float(jn-jic))/float(nn_e)
                     IF (za1 < 1._wp) THEN
                       zwgt1(jn) = 1._wp
@@ -868,5 +868,5 @@
          jpit = jic
       ENDIF
-    
+
       ! Set secondary weights
       DO jn = 1, jpit
@@ -897,24 +897,24 @@
       !!----------------------------------------------------------------------
       !
-      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
+      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
          !                                   ! ---------------
          IF( ln_rstart .AND. ln_bt_fw .AND. (.NOT.l_1st_euler) ) THEN    !* Read the restart file
-            CALL iom_get( numror, jpdom_auto, 'ub2_b'  , ub2_b  (:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )   
-            CALL iom_get( numror, jpdom_auto, 'vb2_b'  , vb2_b  (:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios ) 
-            CALL iom_get( numror, jpdom_auto, 'un_bf'  , un_bf  (:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )   
-            CALL iom_get( numror, jpdom_auto, 'vn_bf'  , vn_bf  (:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios ) 
+            CALL iom_get( numror, jpdom_auto, 'ub2_b'  , ub2_b  (:,:), cd_type = 'U', psgn = -1._wp )
+            CALL iom_get( numror, jpdom_auto, 'vb2_b'  , vb2_b  (:,:), cd_type = 'V', psgn = -1._wp )
+            CALL iom_get( numror, jpdom_auto, 'un_bf'  , un_bf  (:,:), cd_type = 'U', psgn = -1._wp )
+            CALL iom_get( numror, jpdom_auto, 'vn_bf'  , vn_bf  (:,:), cd_type = 'V', psgn = -1._wp )
             IF( .NOT.ln_bt_av ) THEN
-               CALL iom_get( numror, jpdom_auto, 'sshbb_e'  , sshbb_e(:,:), cd_type = 'T', psgn =  1._wp, ldxios = lrxios )   
-               CALL iom_get( numror, jpdom_auto, 'ubb_e'    ,   ubb_e(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )   
-               CALL iom_get( numror, jpdom_auto, 'vbb_e'    ,   vbb_e(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
-               CALL iom_get( numror, jpdom_auto, 'sshb_e'   ,  sshb_e(:,:), cd_type = 'T', psgn =  1._wp, ldxios = lrxios ) 
-               CALL iom_get( numror, jpdom_auto, 'ub_e'     ,    ub_e(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )   
-               CALL iom_get( numror, jpdom_auto, 'vb_e'     ,    vb_e(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
+               CALL iom_get( numror, jpdom_auto, 'sshbb_e'  , sshbb_e(:,:), cd_type = 'T', psgn =  1._wp )
+               CALL iom_get( numror, jpdom_auto, 'ubb_e'    ,   ubb_e(:,:), cd_type = 'U', psgn = -1._wp )
+               CALL iom_get( numror, jpdom_auto, 'vbb_e'    ,   vbb_e(:,:), cd_type = 'V', psgn = -1._wp )
+               CALL iom_get( numror, jpdom_auto, 'sshb_e'   ,  sshb_e(:,:), cd_type = 'T', psgn =  1._wp )
+               CALL iom_get( numror, jpdom_auto, 'ub_e'     ,    ub_e(:,:), cd_type = 'U', psgn = -1._wp )
+               CALL iom_get( numror, jpdom_auto, 'vb_e'     ,    vb_e(:,:), cd_type = 'V', psgn = -1._wp )
             ENDIF
 #if defined key_agrif
             ! Read time integrated fluxes
             IF ( .NOT.Agrif_Root() ) THEN
-               CALL iom_get( numror, jpdom_auto, 'ub2_i_b'  , ub2_i_b(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )   
-               CALL iom_get( numror, jpdom_auto, 'vb2_i_b'  , vb2_i_b(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
+               CALL iom_get( numror, jpdom_auto, 'ub2_i_b'  , ub2_i_b(:,:), cd_type = 'U', psgn = -1._wp )
+               CALL iom_get( numror, jpdom_auto, 'vb2_i_b'  , vb2_i_b(:,:), cd_type = 'V', psgn = -1._wp )
             ELSE
                ub2_i_b(:,:) = 0._wp   ;   vb2_i_b(:,:) = 0._wp   ! used in the 1st update of agrif
@@ -935,26 +935,24 @@
          !                                   ! -------------------
          IF(lwp) WRITE(numout,*) '---- ts_rst ----'
-         IF( lwxios ) CALL iom_swap(      cwxios_context          )
-         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:), ldxios = lwxios )
-         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:), ldxios = lwxios )
-         CALL iom_rstput( kt, nitrst, numrow, 'un_bf'   , un_bf  (:,:), ldxios = lwxios )
-         CALL iom_rstput( kt, nitrst, numrow, 'vn_bf'   , vn_bf  (:,:), ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:) )
+         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:) )
+         CALL iom_rstput( kt, nitrst, numrow, 'un_bf'   , un_bf  (:,:) )
+         CALL iom_rstput( kt, nitrst, numrow, 'vn_bf'   , vn_bf  (:,:) )
          !
          IF (.NOT.ln_bt_av) THEN
-            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:), ldxios = lwxios ) 
-            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:), ldxios = lwxios )
-            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:), ldxios = lwxios )
-            CALL iom_rstput( kt, nitrst, numrow, 'sshb_e'   ,  sshb_e(:,:), ldxios = lwxios )
-            CALL iom_rstput( kt, nitrst, numrow, 'ub_e'     ,    ub_e(:,:), ldxios = lwxios )
-            CALL iom_rstput( kt, nitrst, numrow, 'vb_e'     ,    vb_e(:,:), ldxios = lwxios )
+            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:) )
+            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:) )
+            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:) )
+            CALL iom_rstput( kt, nitrst, numrow, 'sshb_e'   ,  sshb_e(:,:) )
+            CALL iom_rstput( kt, nitrst, numrow, 'ub_e'     ,    ub_e(:,:) )
+            CALL iom_rstput( kt, nitrst, numrow, 'vb_e'     ,    vb_e(:,:) )
          ENDIF
 #if defined key_agrif
          ! Save time integrated fluxes
          IF ( .NOT.Agrif_Root() ) THEN
-            CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b'  , ub2_i_b(:,:), ldxios = lwxios )
-            CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b'  , vb2_i_b(:,:), ldxios = lwxios )
+            CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b'  , ub2_i_b(:,:) )
+            CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b'  , vb2_i_b(:,:) )
          ENDIF
 #endif
-         IF( lwxios ) CALL iom_swap(      cxios_context          )
       ENDIF
       !
@@ -986,5 +984,5 @@
       ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
       IF( ln_bt_auto )   nn_e = CEILING( rn_Dt / rn_bt_cmax * zcmax)
-      
+
       rDt_e = rn_Dt / REAL( nn_e , wp )
       zcmax = zcmax * rDt_e
@@ -1022,5 +1020,5 @@
          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '           Dirac'
          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_e'
-         CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_e' 
+         CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_e'
          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
       END SELECT
@@ -1040,5 +1038,5 @@
       ENDIF
       IF( zcmax>0.9_wp ) THEN
-         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' )          
+         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' )
       ENDIF
       !
@@ -1049,31 +1047,7 @@
       CALL ts_rst( nit000, 'READ' )
       !
-      IF( lwxios ) THEN
-! define variables in restart file when writing with XIOS
-         CALL iom_set_rstw_var_active('ub2_b')
-         CALL iom_set_rstw_var_active('vb2_b')
-         CALL iom_set_rstw_var_active('un_bf')
-         CALL iom_set_rstw_var_active('vn_bf')
-         !
-         IF (.NOT.ln_bt_av) THEN
-            CALL iom_set_rstw_var_active('sshbb_e')
-            CALL iom_set_rstw_var_active('ubb_e')
-            CALL iom_set_rstw_var_active('vbb_e')
-            CALL iom_set_rstw_var_active('sshb_e')
-            CALL iom_set_rstw_var_active('ub_e')
-            CALL iom_set_rstw_var_active('vb_e')
-         ENDIF
-#if defined key_agrif
-         ! Save time integrated fluxes
-         IF ( .NOT.Agrif_Root() ) THEN
-            CALL iom_set_rstw_var_active('ub2_i_b')
-            CALL iom_set_rstw_var_active('vb2_i_b')
-         ENDIF
-#endif
-      ENDIF
-      !
    END SUBROUTINE dyn_spg_ts_init
 
-   
+
    SUBROUTINE dyn_cor_2D_init( Kmm )
       !!---------------------------------------------------------------------
@@ -1102,5 +1076,5 @@
             DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                    &
-                    &           ht(ji  ,jj  ) + ht(ji+1,jj  )   ) * 0.25_wp  
+                    &           ht(ji  ,jj  ) + ht(ji+1,jj  )   ) * 0.25_wp
                IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
             END_2D
@@ -1138,17 +1112,17 @@
          zwz(:,:) = 0._wp
          zhf(:,:) = 0._wp
-         
-         !!gm  assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed 
+
+         !!gm  assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed
 !!gm    A priori a better value should be something like :
-!!gm          zhf(i,j) = masked sum of  ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1) 
-!!gm                     divided by the sum of the corresponding mask 
-!!gm 
-!!            
+!!gm          zhf(i,j) = masked sum of  ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)
+!!gm                     divided by the sum of the corresponding mask
+!!gm
+!!
          IF( .NOT.ln_sco ) THEN
-  
+
    !!gm  agree the JC comment  : this should be done in a much clear way
-  
+
    ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
-   !     Set it to zero for the time being 
+   !     Set it to zero for the time being
    !              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level
    !              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
@@ -1177,5 +1151,5 @@
          END DO
          CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )
-         ! JC: TBC. hf should be greater than 0 
+         ! JC: TBC. hf should be greater than 0
          DO_2D( 1, 1, 1, 1 )
             IF( zhf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zhf(ji,jj)
@@ -1183,5 +1157,5 @@
          zwz(:,:) = ff_f(:,:) * zwz(:,:)
       END SELECT
-      
+
    END SUBROUTINE dyn_cor_2d_init
 
@@ -1209,8 +1183,8 @@
                !
             zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv                    &
-               &               * (  e1e2t(ji,jj+1)*pht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) )   & 
-               &                  + e1e2t(ji,jj  )*pht(ji,jj  )*ff_t(ji,jj  ) * ( punb(ji,jj  ) + punb(ji-1,jj  ) )   ) 
-         END_2D
-         !         
+               &               * (  e1e2t(ji,jj+1)*pht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) )   &
+               &                  + e1e2t(ji,jj  )*pht(ji,jj  )*ff_t(ji,jj  ) * ( punb(ji,jj  ) + punb(ji-1,jj  ) )   )
+         END_2D
+         !
       CASE( np_ENE , np_MIX )        ! energy conserving scheme (t-point) ENE or MIX
          DO_2D( 0, 0, 0, 0 )
@@ -1234,5 +1208,5 @@
          END_2D
          !
-      CASE( np_EET , np_EEN )      ! energy & enstrophy scheme (using e3t or e3f)         
+      CASE( np_EET , np_EEN )      ! energy & enstrophy scheme (using e3t or e3f)
          DO_2D( 0, 0, 0, 0 )
             zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zhV(ji  ,jj  ) &
@@ -1254,9 +1228,9 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE wad_lmt  ***
-      !!                    
-      !! ** Purpose :   set wetting & drying mask at tracer points 
-      !!              for the current barotropic sub-step 
-      !!
-      !! ** Method  :   ??? 
+      !!
+      !! ** Purpose :   set wetting & drying mask at tracer points
+      !!              for the current barotropic sub-step
+      !!
+      !! ** Method  :   ???
       !!
       !! ** Action  :  ptmsk : wetting & drying t-mask
@@ -1268,16 +1242,16 @@
       !!----------------------------------------------------------------------
       !
-      IF( ln_wd_dl_rmp ) THEN     
+      IF( ln_wd_dl_rmp ) THEN
          DO_2D( 1, 1, 1, 1 )
-            IF    ( pssh(ji,jj) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN 
-               !           IF    ( pssh(ji,jj) + ht_0(ji,jj) >          rn_wdmin2 ) THEN 
+            IF    ( pssh(ji,jj) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN
+               !           IF    ( pssh(ji,jj) + ht_0(ji,jj) >          rn_wdmin2 ) THEN
                ptmsk(ji,jj) = 1._wp
             ELSEIF( pssh(ji,jj) + ht_0(ji,jj) >          rn_wdmin1 ) THEN
                ptmsk(ji,jj) = TANH( 50._wp*( ( pssh(ji,jj) + ht_0(ji,jj) -  rn_wdmin1 )*r_rn_wdmin1) )
-            ELSE 
+            ELSE
                ptmsk(ji,jj) = 0._wp
             ENDIF
          END_2D
-      ELSE  
+      ELSE
          DO_2D( 1, 1, 1, 1 )
             IF ( pssh(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN   ;   ptmsk(ji,jj) = 1._wp
@@ -1293,9 +1267,9 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE wad_lmt  ***
-      !!                    
-      !! ** Purpose :   set wetting & drying mask at tracer points 
-      !!              for the current barotropic sub-step 
-      !!
-      !! ** Method  :   ??? 
+      !!
+      !! ** Purpose :   set wetting & drying mask at tracer points
+      !!              for the current barotropic sub-step
+      !!
+      !! ** Method  :   ???
       !!
       !! ** Action  :  ptmsk : wetting & drying t-mask
@@ -1309,6 +1283,6 @@
       !
       DO_2D( 1, 1, 1, 0 )   ! not jpi-column
-         IF ( phU(ji,jj) > 0._wp ) THEN   ;   pUmsk(ji,jj) = pTmsk(ji  ,jj) 
-         ELSE                             ;   pUmsk(ji,jj) = pTmsk(ji+1,jj)  
+         IF ( phU(ji,jj) > 0._wp ) THEN   ;   pUmsk(ji,jj) = pTmsk(ji  ,jj)
+         ELSE                             ;   pUmsk(ji,jj) = pTmsk(ji+1,jj)
          ENDIF
          phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj)
@@ -1318,7 +1292,7 @@
       DO_2D( 1, 0, 1, 1 )   ! not jpj-row
          IF ( phV(ji,jj) > 0._wp ) THEN   ;   pVmsk(ji,jj) = pTmsk(ji,jj  )
-         ELSE                             ;   pVmsk(ji,jj) = pTmsk(ji,jj+1)  
-         ENDIF
-         phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj) 
+         ELSE                             ;   pVmsk(ji,jj) = pTmsk(ji,jj+1)
+         ENDIF
+         phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj)
          pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj)
       END_2D
@@ -1331,5 +1305,5 @@
       !!                   ***  ROUTINE  wad_sp  ***
       !!
-      !! ** Purpose : 
+      !! ** Purpose :
       !!----------------------------------------------------------------------
       INTEGER  ::   ji ,jj               ! dummy loop indices
@@ -1364,5 +1338,5 @@
               &      MAX(   pshn(ji,jj)              ,  pshn(ji,jj+1) ) >                &
               &      MAX(  -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
-         
+
          IF(ll_tmp1) THEN
             zcpy(ji,jj) = 1.0_wp
@@ -1376,7 +1350,7 @@
          ENDIF
       END_2D
-            
+
    END SUBROUTINE wad_spg
-     
+
 
 
@@ -1384,10 +1358,10 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_drg_init  ***
-      !!                    
-      !! ** Purpose : - add the baroclinic top/bottom drag contribution to 
+      !!
+      !! ** Purpose : - add the baroclinic top/bottom drag contribution to
       !!              the baroclinic part of the barotropic RHS
       !!              - compute the barotropic drag coefficients
       !!
-      !! ** Method  :   computation done over the INNER domain only 
+      !! ** Method  :   computation done over the INNER domain only
       !!----------------------------------------------------------------------
       INTEGER                             , INTENT(in   ) ::  Kbb, Kmm           ! ocean time level indices
@@ -1406,5 +1380,5 @@
       !
       IF( ln_isfcav.OR.ln_drgice_imp ) THEN          ! top+bottom friction (ocean cavities)
-         
+
          DO_2D( 0, 0, 0, 0 )
             pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
@@ -1421,16 +1395,16 @@
       !
       IF( ln_bt_fw ) THEN                 ! FORWARD integration: use NOW bottom baroclinic velocities
-         
+
          DO_2D( 0, 0, 0, 0 )
-            ikbu = mbku(ji,jj)       
-            ikbv = mbkv(ji,jj)    
+            ikbu = mbku(ji,jj)
+            ikbv = mbkv(ji,jj)
             zu_i(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm)
             zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm)
          END_2D
       ELSE                                ! CENTRED integration: use BEFORE bottom baroclinic velocities
-         
+
          DO_2D( 0, 0, 0, 0 )
-            ikbu = mbku(ji,jj)       
-            ikbv = mbkv(ji,jj)    
+            ikbu = mbku(ji,jj)
+            ikbv = mbkv(ji,jj)
             zu_i(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb)
             zv_i(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb)
@@ -1441,11 +1415,11 @@
          zztmp = -1._wp / rDt_e
          DO_2D( 0, 0, 0, 0 )
-            pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) *  wdrampu(ji,jj) * MAX(                                 & 
+            pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) *  wdrampu(ji,jj) * MAX(                                 &
                  &                              r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp  )
-            pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) *  wdrampv(ji,jj) * MAX(                                 & 
+            pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) *  wdrampv(ji,jj) * MAX(                                 &
                  &                              r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp  )
          END_2D
       ELSE                    ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
-         
+
          DO_2D( 0, 0, 0, 0 )
             pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj)
@@ -1459,5 +1433,5 @@
          !
          IF( ln_bt_fw ) THEN                ! FORWARD integration: use NOW top baroclinic velocity
-            
+
             DO_2D( 0, 0, 0, 0 )
                iktu = miku(ji,jj)
@@ -1467,5 +1441,5 @@
             END_2D
          ELSE                                ! CENTRED integration: use BEFORE top baroclinic velocity
-            
+
             DO_2D( 0, 0, 0, 0 )
                iktu = miku(ji,jj)
@@ -1477,5 +1451,5 @@
          !
          !                    ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
-         
+
          DO_2D( 0, 0, 0, 0 )
             pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj)
@@ -1498,6 +1472,6 @@
       !                             ! set Half-step back interpolation coefficient
       IF    ( jn==1 .AND. ll_init ) THEN   !* Forward-backward
-         za0 = 1._wp                        
-         za1 = 0._wp                           
+         za0 = 1._wp
+         za1 = 0._wp
          za2 = 0._wp
          za3 = 0._wp
@@ -1506,7 +1480,7 @@
          za1 =-0.1666666666666_wp                 ! za1 = gam
          za2 = 0.0833333333333_wp                 ! za2 = eps
-         za3 = 0._wp              
-      ELSE                                 !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 
-         IF( rn_bt_alpha == 0._wp ) THEN      ! Time diffusion  
+         za3 = 0._wp
+      ELSE                                 !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
+         IF( rn_bt_alpha == 0._wp ) THEN      ! Time diffusion
             za0 = 0.614_wp                        ! za0 = 1/2 +   gam + 2*eps
             za1 = 0.285_wp                        ! za1 = 1/2 - 2*gam - 3*eps
@@ -1520,5 +1494,5 @@
             za2 = zgamma
             za3 = zepsilon
-         ENDIF 
+         ENDIF
       ENDIF
    END SUBROUTINE ts_bck_interp
