New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 4241 – NEMO

Changeset 4241


Ignore:
Timestamp:
2013-11-19T09:04:21+01:00 (10 years ago)
Author:
acc
Message:

Branch 2013/dev_r3858_NOC_ZTC, #863. Port of Jeromes version of dynspg_ts.F90 into z-tilde framework

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r4194 r4241  
    99   !!             3.3  ! 2010-09  (D. Storkey, E. O'Dea) update for BDY for Shelf configurations 
    1010   !!             3.3  ! 2011-03  (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b 
     11   !!             3.5  ! 2013-07  (J. Chanut) Switch to Forward-backward time stepping 
     12   !!             3.6  ! 2013-11  (A. Coward) Update for z-tilde compatibility 
    1113   !!--------------------------------------------------------------------- 
    1214#if defined key_dynspg_ts   ||   defined key_esopa 
     
    1618   !!   dyn_spg_ts  : compute surface pressure gradient trend using a time- 
    1719   !!                 splitting scheme and add to the general trend  
    18    !!   ts_rst      : read/write the time-splitting restart fields in the ocean restart file 
    1920   !!---------------------------------------------------------------------- 
    2021   USE oce             ! ocean dynamics and tracers 
     
    2425   USE phycst          ! physical constants 
    2526   USE domvvl          ! variable volume 
    26    USE zdfbfr          ! bottom friction 
    2727   USE dynvor          ! vorticity term 
    28    USE dynadv          ! advection 
    29    USE obc_oce         ! Lateral open boundary condition 
    30    USE obc_par         ! open boundary condition parameters 
    31    USE obcdta          ! open boundary condition data      
    32    USE obcfla          ! Flather open boundary condition   
    3328   USE bdy_par         ! for lk_bdy 
    3429   USE bdy_oce         ! Lateral open boundary condition 
    35    USE bdydta          ! open boundary condition data      
     30   USE bdytides        ! open boundary condition data      
    3631   USE bdydyn2d        ! open boundary conditions on barotropic variables 
    37    USE bdytides        !  
    38    USE sbctide 
    39    USE updtide 
     32   USE sbctide         ! tides 
     33   USE updtide         ! tide potential 
    4034   USE lib_mpp         ! distributed memory computing library 
    4135   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     
    4337   USE in_out_manager  ! I/O manager 
    4438   USE iom             ! IOM library 
     39   USE restart         ! only for lrst_oce 
    4540   USE zdf_oce         ! Vertical diffusion 
    4641   USE wrk_nemo        ! Memory Allocation 
    47    USE timing          ! Timing 
     42   USE timing          ! Timing     
     43   USE sbcapr          ! surface boundary condition: atmospheric pressure 
     44   USE dynadv, ONLY: ln_dynadv_vec 
    4845 
    4946 
    5047   IMPLICIT NONE 
    5148   PRIVATE 
     49 
     50   PUBLIC dyn_spg_ts        ! routine called in dynspg.F90  
     51   PUBLIC dyn_spg_ts_alloc  !    "      "     "    " 
     52   PUBLIC dyn_spg_ts_init   !    "      "     "    " 
    5253 
    5354   ! Potential namelist parameters below to be read in dyn_spg_ts_init 
     
    6566                    wgtbtp1, &              ! Primary weights used for time filtering of barotropic variables 
    6667                    wgtbtp2                 ! Secondary weights used for time filtering of barotropic variables 
    67    ! 
    68    PUBLIC dyn_spg_ts        ! routine called by step.F90 
    69    PUBLIC ts_rst            ! routine called by istate.F90 
    70    PUBLIC dyn_spg_ts_alloc  ! routine called by dynspg.F90 
    71    PUBLIC dyn_spg_ts_init   !    "      "     "    " 
    72  
    7368 
    7469   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zwz          ! ff/h at F points 
     
    9691#  include "vectopt_loop_substitute.h90" 
    9792   !!---------------------------------------------------------------------- 
    98    !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
    99    !! $Id$ 
     93   !! NEMO/OPA 3.5 , NEMO Consortium (2013) 
     94   !! $Id: dynspg_ts.F90 
    10095   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    10196   !!---------------------------------------------------------------------- 
     
    132127   END FUNCTION dyn_spg_ts_alloc 
    133128 
    134  
    135129   SUBROUTINE dyn_spg_ts( kt ) 
    136130      !!---------------------------------------------------------------------- 
    137       !!                  ***  routine dyn_spg_ts  *** 
    138131      !! 
    139       !! ** Purpose :   Compute the now trend due to the surface pressure 
    140       !!      gradient in case of free surface formulation with time-splitting. 
    141       !!      Add it to the general trend of momentum equation. 
     132      !! ** Purpose :    
     133      !!      -Compute the now trend due to the explicit time stepping 
     134      !!      of the quasi-linear barotropic system. Barotropic variables are 
     135      !!      advanced from internal time steps "n" to "n+1" (if ln_bt_cen=F) 
     136      !!      or from "n-1" to "n+1" time steps (if ln_bt_cen=T) with a 
     137      !!      generalized forward-backward (see ref. below) time stepping. 
     138      !!      -Update the free surface at step "n+1" (ssha, zsshu_a, zsshv_a). 
     139      !!      -Compute barotropic advective velocities at step "n" to be used  
     140      !!      to advect tracers latter on. These are compliant with discrete 
     141      !!      continuity equation taken at the baroclinic time steps, thus  
     142      !!      ensuring tracers conservation. 
    142143      !! 
    143       !! ** Method  :   Free surface formulation with time-splitting 
    144       !!      -1- Save the vertically integrated trend. This general trend is 
    145       !!          held constant over the barotropic integration. 
    146       !!          The Coriolis force is removed from the general trend as the 
    147       !!          surface gradient and the Coriolis force are updated within 
    148       !!          the barotropic integration. 
    149       !!      -2- Barotropic loop : updates of sea surface height (ssha_e) and  
    150       !!          barotropic velocity (ua_e and va_e) through barotropic  
    151       !!          momentum and continuity integration. Barotropic former  
    152       !!          variables are time averaging over the full barotropic cycle 
    153       !!          (= 2 * baroclinic time step) and saved in uX_b  
    154       !!          and vX_b (X specifying after, now or before). 
    155       !!      -3- The new general trend becomes : 
    156       !!          ua = ua - sum_k(ua)/H + ( un_b - ub_b ) 
     144      !! ** Method  :   
    157145      !! 
    158       !! ** Action : - Update (ua,va) with the surf. pressure gradient trend 
     146      !! ** Action : - Update barotropic velocities: ua_b, va_b 
     147      !!             - Update trend (ua,va) with barotropic component 
     148      !!             - Update ssha, zsshu_a, zsshv_a 
     149      !!             - Update barotropic advective velocity at kt=now 
    159150      !! 
    160       !! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL 
     151      !! References : Shchepetkin, A.F. and J.C. McWilliams, 2005:  
     152      !!              The regional oceanic modeling system (ROMS):  
     153      !!              a split-explicit, free-surface, 
     154      !!              topography-following-coordinate oceanic model.  
     155      !!              Ocean Modelling, 9, 347-404.  
    161156      !!--------------------------------------------------------------------- 
    162157      ! 
    163158      INTEGER, INTENT(in)  ::   kt   ! ocean time-step index 
    164159      ! 
    165       LOGICAL  ::   ll_fw_start      ! if true, forward integration  
    166       LOGICAL  ::   ll_init          ! if true, special startup of 2d equations 
    167       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    168 !     INTEGER  ::   icycle           ! local scalar 
    169       INTEGER  ::   ioffset          ! local scalar 
    170       INTEGER  ::   ikbu, ikbv       ! local scalar 
    171       REAL(wp) ::   zraur, zcoef, z2dt_e, z1_2dt_b, z2dt_bf   ! local scalars 
    172       REAL(wp) ::   zx1, zy1, zx2, zy2                        !   -      - 
    173       REAL(wp) ::   z1_12, z1_8, z1_4, z1_2 
    174       REAL(wp) ::   za0, za1, za2, za3                        !   -      - 
    175       REAL(wp) ::   zu_spg, zu_cor, zu_sld, zu_asp            !   -      - 
    176       REAL(wp) ::   zv_spg, zv_cor, zv_sld, zv_asp            !   -      - 
    177       REAL(wp) ::   ua_btm, va_btm                            !   -      - 
    178       ! 
    179       REAL(wp), POINTER, DIMENSION(:,:) :: zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv  
    180       REAL(wp), POINTER, DIMENSION(:,:) :: zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e  
    181       REAL(wp), POINTER, DIMENSION(:,:) :: zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum 
    182       REAL(wp), POINTER, DIMENSION(:,:) :: zhu_b, zhv_b 
    183       REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zsshp2_e 
    184       REAL(wp), POINTER, DIMENSION(:,:) :: zhust_e, zhvst_e 
     160      LOGICAL  ::   ll_fw_start        ! if true, forward integration  
     161      LOGICAL  ::   ll_init         ! if true, special startup of 2d equations 
     162      INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices 
     163      INTEGER  ::   ikbu, ikbv, noffset      ! local integers 
     164      REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf    ! local scalars 
     165      REAL(wp) ::   zx1, zy1, zx2, zy2          !   -      - 
     166      REAL(wp) ::   z1_12, z1_8, z1_4, z1_2  !   -      - 
     167      REAL(wp) ::   zu_spg, zv_spg     !   -      - 
     168      REAL(wp) ::   zhura, zhvra          !   -      - 
     169      REAL(wp) ::   za0, za1, za2, za3    !   -      - 
     170      ! 
     171      REAL(wp), POINTER, DIMENSION(:,:) :: zun_e, zvn_e, zsshp2_e 
     172      REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc 
     173      REAL(wp), POINTER, DIMENSION(:,:) :: zu_sum, zv_sum, zwx, zwy, zhdiv 
     174      REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e 
     175      REAL(wp), POINTER, DIMENSION(:,:) :: zhur_b, zhvr_b 
     176      REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 
     177      REAL(wp), POINTER, DIMENSION(:,:) :: zht, zhf 
    185178      !!---------------------------------------------------------------------- 
    186179      ! 
    187180      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_ts') 
    188181      ! 
    189       CALL wrk_alloc( jpi, jpj, zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv     ) 
    190       CALL wrk_alloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e   ) 
    191       CALL wrk_alloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum ) 
    192       CALL wrk_alloc( jpi, jpj, zhu_b, zhv_b                                     ) 
    193       CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zsshp2_e                       ) 
    194       CALL wrk_alloc( jpi, jpj, zhust_e, zhvst_e                                 ) 
    195       ! 
    196  
    197       !                                                     !* Local constant initialization 
    198       z1_12    = 1._wp / 12._wp 
    199       z1_8     = 0.125_wp                                   ! coefficients for vorticity estimates 
    200       z1_4     = 0.25_wp         
    201       z1_2     = 0.5_wp         
    202       zraur    = 1._wp / rau0                               ! 1 / volumetric mass 
    203       ! 
    204       zhdiv(:,:) = 0._wp                                    ! barotropic divergence 
    205       zu_sld = 0._wp   ;   zu_asp = 0._wp                   ! tides trends (lk_tide=F) 
    206       zv_sld = 0._wp   ;   zv_asp = 0._wp 
    207  
    208       IF( kt == nit000 .AND. neuler == 0) THEN              ! for implicit bottom friction 
    209         z2dt_bf = rdt                                       ! baroclinic time step (euler timestep) 
    210       ELSE 
    211         z2dt_bf = 2.0_wp * rdt                              ! baroclinic time step 
    212       ENDIF 
    213       z1_2dt_b = 1.0_wp / z2dt_bf                           ! reciprocal of baroclinic time step 
    214       ! 
    215       ll_init = ln_bt_av                                    ! if no time averaging, then no specific restart  
     182      !                                         !* Allocate temporay arrays 
     183      CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv ) 
     184      CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e  ) 
     185      CALL wrk_alloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc) 
     186      CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 
     187      CALL wrk_alloc( jpi, jpj, zhur_b, zhvr_b                                     ) 
     188      CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a                                   ) 
     189      CALL wrk_alloc( jpi, jpj, zht, zhf ) 
     190      ! 
     191      !                                         !* Local constant initialization 
     192      z1_12 = 1._wp / 12._wp  
     193      z1_8  = 0.125_wp                                    
     194      z1_4  = 0.25_wp 
     195      z1_2  = 0.5_wp      
     196      zraur = 1._wp / rau0 
     197      ! 
     198      IF( kt == nit000 .AND. neuler == 0 ) THEN    ! reciprocal of baroclinic time step  
     199        z2dt_bf = rdt 
     200      ELSE 
     201        z2dt_bf = 2.0_wp * rdt 
     202      ENDIF 
     203      z1_2dt_b = 1.0_wp / z2dt_bf  
     204      ! 
     205      ll_init = ln_bt_av                           ! if no time averaging, then no specific restart  
    216206      ll_fw_start = .FALSE. 
    217207      ! 
    218       ! time offset in steps for bdy data update 
    219       IF (.NOT.ln_bt_fw) THEN ; ioffset=-2*nn_baro ; ELSE ;  ioffset = 0 ; ENDIF 
    220       ! 
    221       IF( kt == nit000 ) THEN                   !* initialisation 
     208                                                       ! time offset in steps for bdy data update 
     209      IF (.NOT.ln_bt_fw) THEN ; noffset=-2*nn_baro ; ELSE ;  noffset = 0 ; ENDIF 
     210      ! 
     211      IF( kt == nit000 ) THEN                !* initialisation 
    222212         ! 
    223213         IF(lwp) WRITE(numout,*) 
     
    230220         IF (ln_bt_fw.OR.(neuler==0)) THEN 
    231221           ll_fw_start=.TRUE. 
    232            ioffset = 0 
     222           noffset = 0 
    233223         ELSE 
    234224           ll_fw_start=.FALSE. 
     
    238228         CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2) 
    239229         ! 
    240          IF ((neuler/=0).AND.(ln_bt_fw)) CALL ts_rst( nit000, 'READ' )   ! read or initialize the following fields: un_b, vn_b 
    241          ! 
    242          ua_e  (:,:) = un_b (:,:) 
    243          va_e  (:,:) = vn_b (:,:) 
    244          hu_e  (:,:) = hu   (:,:) 
    245          hv_e  (:,:) = hv   (:,:) 
    246          hur_e (:,:) = hur  (:,:) 
    247          hvr_e (:,:) = hvr  (:,:) 
    248          IF( ln_dynvor_een ) THEN 
    249             ftne(1,:) = 0._wp   ;   ftnw(1,:) = 0._wp   ;   ftse(1,:) = 0._wp   ;   ftsw(1,:) = 0._wp 
     230         IF ((neuler/=0).AND.(ln_bt_fw)) CALL ts_rst( nit000, 'READ' )  
     231         ! 
     232      ENDIF 
     233      ! 
     234      ! Set arrays to remove/compute coriolis trend. 
     235      ! Do it once at kt=nit000 if volume is fixed, else at each long time step. 
     236      ! Note that these arrays are also used during barotropic loop. These are however frozen 
     237      ! although they should be updated in variable volume case. Not a big approximation. 
     238      ! To remove this approximation, copy lines below inside barotropic loop 
     239      ! and update depths at T-F points (ht and hf resp.) at each barotropic time step 
     240      ! 
     241      IF ( kt == nit000 .OR. lk_vvl ) THEN 
     242         IF ( ln_dynvor_een ) THEN 
     243            zht(:,:) = 0. 
     244            DO jk = 1, jpk 
     245               zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 
     246            END DO 
     247            DO jj = 1, jpjm1 
     248               DO ji = 1, jpim1 
     249                  zwz(ji,jj) =   ( zht(ji  ,jj+1) + zht(ji+1,jj+1) +                     & 
     250                        &          zht(ji  ,jj  ) + zht(ji+1,jj  )   )                   & 
     251                        &      / ( MAX( 1.0_wp, tmask(ji  ,jj+1, 1) + tmask(ji+1,jj+1, 1) +    & 
     252                        &                       tmask(ji  ,jj  , 1) + tmask(ji+1,jj  , 1) ) ) 
     253                  IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zwz(ji,jj) 
     254               END DO 
     255            END DO 
     256            CALL lbc_lnk( zwz, 'F', 1._wp ) 
     257            zwz(:,:) = ff(:,:) * zwz(:,:) 
     258 
     259            ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 
    250260            DO jj = 2, jpj 
    251261               DO ji = fs_2, jpi   ! vector opt. 
    252                   ftne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3._wp 
    253                   ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3._wp 
    254                   ftse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3._wp 
    255                   ftsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3._wp 
    256                END DO 
    257             END DO 
    258          ENDIF 
    259          ! 
     262                  ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1) 
     263                  ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  ) 
     264                  ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1) 
     265                  ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  ) 
     266               END DO 
     267            END DO 
     268         ELSE 
     269            zwz(:,:) = 0._wp 
     270            zht(:,:) = 0. 
     271            IF ( .not. ln_sco ) THEN 
     272!              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level 
     273!              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth 
     274!              ENDIF 
     275!              zht(:,:) = gdepw_0(:,:,jk+1) 
     276            ELSE 
     277               zht(:,:) = hbatf(:,:) 
     278            END IF 
     279 
     280            DO jj = 1, jpjm1 
     281               zht(:,jj) = zht(:,jj)*(1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 
     282            END DO 
     283 
     284            DO jk = 1, jpkm1 
     285               DO jj = 1, jpjm1 
     286                  zht(:,jj) = zht(:,jj) + fse3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 
     287               END DO 
     288            END DO 
     289            CALL lbc_lnk( zht, 'F', 1._wp ) 
     290            ! JC: TBC. hf should be greater than 0  
     291            DO jj = 1, jpj 
     292               DO ji = 1, jpi 
     293                  IF( zht(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zht(ji,jj) ! zht is actually hf here but it saves an array 
     294               END DO 
     295            END DO 
     296            zwz(:,:) = ff(:,:) * zwz(:,:) 
     297         ENDIF 
    260298      ENDIF 
    261299      ! 
     
    267305      ENDIF 
    268306 
     307      ! before inverse water column height at u- and v- points 
     308      IF( lk_vvl ) THEN 
     309         zhur_b(:,:) = 0. 
     310         zhvr_b(:,:) = 0. 
     311         DO jk = 1, jpk 
     312            zhur_b(:,:) = zhur_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk) 
     313            zhvr_b(:,:) = zhvr_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 
     314         END DO 
     315         zhur_b(:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1. - umask(:,:,1) ) 
     316         zhvr_b(:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1. - vmask(:,:,1) ) 
     317      ELSE 
     318         zhur_b(:,:) = hur(:,:) 
     319         zhvr_b(:,:) = hvr(:,:) 
     320      ENDIF 
     321                           
    269322      ! ----------------------------------------------------------------------------- 
    270323      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step) 
    271324      ! ----------------------------------------------------------------------------- 
    272325      !       
     326      ! Some vertical sums (at now and before time steps) below could be suppressed  
     327      ! if one swap barotropic arrays somewhere 
     328      ! 
    273329      !                                   !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated) 
    274       !                                   ! -------------------------- 
    275       zua(:,:) = 0._wp   ;   zun(:,:) = 0._wp   ;   ub_b(:,:) = 0._wp 
    276       zva(:,:) = 0._wp   ;   zvn(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp 
     330      !                                   ! -------------------------------------------------- 
     331      zu_frc(:,:) = 0._wp   ;   ub_b(:,:) = 0._wp  ;  un_b(:,:) = 0._wp 
     332      zv_frc(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp  ;  vn_b(:,:) = 0._wp 
    277333      ! 
    278334      DO jk = 1, jpkm1 
     
    284340            DO ji = 1, jpi 
    285341#endif 
    286                !                                                                              ! now trend 
    287                zua(ji,jj) = zua(ji,jj) + fse3u_n(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 
    288                zva(ji,jj) = zva(ji,jj) + fse3v_n(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 
    289                !                                                                              ! now velocity  
    290                zun(ji,jj) = zun(ji,jj) + fse3u_n(ji,jj,jk) * un(ji,jj,jk) 
    291                zvn(ji,jj) = zvn(ji,jj) + fse3v_n(ji,jj,jk) * vn(ji,jj,jk)                
    292                ! 
     342               !        ! now trend:                                                                    
     343               zu_frc(ji,jj) = zu_frc(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 
     344               zv_frc(ji,jj) = zv_frc(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)          
     345               !        ! now bt transp:                   
     346               un_b(ji,jj) = un_b(ji,jj) + fse3u(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)         
     347               vn_b(ji,jj) = vn_b(ji,jj) + fse3v(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) 
     348               !  ! before bt transp: 
    293349               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk) 
    294350               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk)  * vmask(ji,jj,jk) 
     
    296352         END DO 
    297353      END DO 
    298  
    299       ! before inverse water column height at u- and v- points 
     354      ! 
     355      zu_frc(:,:) = zu_frc(:,:) * hur(:,:) 
     356      zv_frc(:,:) = zv_frc(:,:) * hvr(:,:) 
     357      ! 
    300358      IF( lk_vvl ) THEN 
    301          zhu_b(:,:) = 0. 
    302          zhv_b(:,:) = 0. 
    303          DO jk = 1, jpk 
    304             zhu_b(:,:) = zhu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk) 
    305             zhv_b(:,:) = zhv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 
    306          END DO 
    307          zhu_b(:,:) = umask(:,:,1) / ( zhu_b(:,:) + 1. - umask(:,:,1) ) 
    308          zhv_b(:,:) = vmask(:,:,1) / ( zhv_b(:,:) + 1. - vmask(:,:,1) ) 
    309       ELSE 
    310          zhu_b(:,:) = hur(:,:) 
    311          zhv_b(:,:) = hvr(:,:) 
    312       ENDIF 
    313  
     359          ub_b(:,:) = ub_b(:,:) * zhur_b(:,:) 
     360          vb_b(:,:) = vb_b(:,:) * zhvr_b(:,:) 
     361      ELSE 
     362          ub_b(:,:) = ub_b(:,:) * hur(:,:) 
     363          vb_b(:,:) = vb_b(:,:) * hvr(:,:) 
     364      ENDIF 
     365      ! 
    314366      !                                   !* baroclinic momentum trend (remove the vertical mean trend) 
    315       DO jk = 1, jpkm1                    ! -------------------------- 
     367      DO jk = 1, jpkm1                    ! ----------------------------------------------------------- 
    316368         DO jj = 2, jpjm1 
    317369            DO ji = fs_2, fs_jpim1   ! vector opt. 
    318                ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj) 
    319                va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj) 
     370               ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk) 
     371               va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk) 
    320372            END DO 
    321373         END DO 
    322374      END DO 
    323  
    324       !                                   !* barotropic Coriolis trends * H (vorticity scheme dependent) 
    325       !                                   ! ---------------------------==== 
    326       zwx(:,:) = zun(:,:) * e2u(:,:)                   ! now transport  
    327       zwy(:,:) = zvn(:,:) * e1v(:,:) 
     375      !                                   !* barotropic Coriolis trends (vorticity scheme dependent) 
     376      !                                   ! -------------------------------------------------------- 
     377      zwx(:,:) = un_b(:,:) * e2u(:,:)           ! now transport  
     378      zwy(:,:) = vn_b(:,:) * e1v(:,:) 
    328379      ! 
    329380      IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme 
     
    335386               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
    336387               ! energy conserving formulation for planetary vorticity term 
    337                zcu(ji,jj) = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 ) 
    338                zcv(ji,jj) =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) 
    339             END DO 
    340          END DO 
    341          ! 
    342       ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme 
     388               zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
     389               zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
     390            END DO 
     391         END DO 
     392         ! 
     393      ELSEIF ( ln_dynvor_ens ) THEN             ! enstrophy conserving scheme 
    343394         DO jj = 2, jpjm1 
    344395            DO ji = fs_2, fs_jpim1   ! vector opt. 
    345                zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
    346                zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
    347                zcu(ji,jj)  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) ) 
    348                zcv(ji,jj)  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) 
    349             END DO 
    350          END DO 
    351          ! 
    352       ELSEIF ( ln_dynvor_een ) THEN                    ! enstrophy and energy conserving scheme 
     396               zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) & 
     397                 &            + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
     398               zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) & 
     399                 &            + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
     400               zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) ) 
     401               zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) ) 
     402            END DO 
     403         END DO 
     404         ! 
     405      ELSEIF ( ln_dynvor_een ) THEN             ! enstrophy and energy conserving scheme 
    353406         DO jj = 2, jpjm1 
    354407            DO ji = fs_2, fs_jpim1   ! vector opt. 
    355                zcu(ji,jj) = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   & 
    356                   &                                + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 
    357                zcv(ji,jj) = - z1_4 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   & 
    358                   &                                + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) 
    359             END DO 
    360          END DO 
    361          ! 
    362       ENDIF 
    363  
     408               zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) & 
     409                &                                      + ftnw(ji+1,jj) * zwy(ji+1,jj  ) & 
     410                &                                      + ftse(ji,jj  ) * zwy(ji  ,jj-1) & 
     411                &                                      + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 
     412               zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 
     413                &                                      + ftse(ji,jj+1) * zwx(ji  ,jj+1) & 
     414                &                                      + ftnw(ji,jj  ) * zwx(ji-1,jj  ) & 
     415                &                                      + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) 
     416            END DO 
     417         END DO 
     418         ! 
     419      ENDIF  
     420      ! 
     421      un_b (:,:) = un_b(:,:) * hur(:,:)         ! Revert now transport to barotropic velocities 
     422      vn_b (:,:) = vn_b(:,:) * hvr(:,:)   
    364423      !                                   !* Right-Hand-Side of the barotropic momentum equation 
    365424      !                                   ! ---------------------------------------------------- 
    366       IF( lk_vvl ) THEN                         ! Variable volume : remove both Coriolis and Surface pressure gradient 
     425      IF( lk_vvl ) THEN                         ! Variable volume : remove surface pressure gradient 
    367426         DO jj = 2, jpjm1  
    368427            DO ji = fs_2, fs_jpim1   ! vector opt. 
    369                zcu(ji,jj) = zcu(ji,jj) - grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn(ji+1,jj  )    & 
    370                   &                              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  )  ) * hu(ji,jj) / e1u(ji,jj) 
    371                zcv(ji,jj) = zcv(ji,jj) - grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn(ji  ,jj+1)    & 
    372                   &                              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  )  ) * hv(ji,jj) / e2v(ji,jj) 
    373             END DO 
    374          END DO 
    375       ENDIF 
    376  
    377       DO jj = 2, jpjm1                             ! Remove coriolis term (and possibly spg) from barotropic trend 
     428               zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) / e1u(ji,jj) 
     429               zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) / e2v(ji,jj) 
     430            END DO 
     431         END DO 
     432      ENDIF 
     433 
     434      DO jj = 2, jpjm1                          ! Remove coriolis term (and possibly spg) from barotropic trend 
    378435         DO ji = fs_2, fs_jpim1 
    379              zua(ji,jj) = zua(ji,jj) - zcu(ji,jj) 
    380              zva(ji,jj) = zva(ji,jj) - zcv(ji,jj) 
     436             zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1) 
     437             zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1) 
    381438          END DO 
    382       END DO 
    383  
    384                      
    385       !                                             ! Remove barotropic contribution of bottom friction  
    386       !                                             ! from the barotropic transport trend 
    387       zcoef = -1._wp * z1_2dt_b 
    388  
    389       IF(ln_bfrimp) THEN 
    390       !                                   ! Remove the bottom stress trend from 3-D sea surface level gradient 
    391       !                                   ! and Coriolis forcing in case of 3D semi-implicit bottom friction  
    392         DO jj = 2, jpjm1          
    393            DO ji = fs_2, fs_jpim1 
    394               ikbu = mbku(ji,jj) 
    395               ikbv = mbkv(ji,jj) 
    396               ua_btm = zcu(ji,jj) * z2dt_bf * hur(ji,jj) * umask (ji,jj,ikbu) 
    397               va_btm = zcv(ji,jj) * z2dt_bf * hvr(ji,jj) * vmask (ji,jj,ikbv) 
    398  
    399               zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * ua_btm 
    400               zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * va_btm 
    401            END DO 
    402         END DO 
    403  
    404       ELSE 
    405  
    406 # if defined key_vectopt_loop 
    407         DO jj = 1, 1 
    408            DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling) 
    409 # else 
    410         DO jj = 2, jpjm1 
    411            DO ji = 2, jpim1 
    412 # endif 
    413             ! Apply stability criteria for bottom friction 
    414             !RBbug for vvl and external mode we may need to use varying fse3 
    415             !!gm  Rq: the bottom e3 present the smallest variation, the use of e3u_0 is not a big approx. 
    416               zbfru(ji,jj) = MAX(  bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef  ) 
    417               zbfrv(ji,jj) = MAX(  bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef  ) 
    418            END DO 
    419         END DO 
    420  
    421         DO jj = 2, jpjm1 
    422            DO ji = fs_2, fs_jpim1   ! vector opt. 
    423               zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * zhu_b(ji,jj) 
    424               zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * zhv_b(ji,jj) 
    425            END DO 
    426         END DO 
    427       END IF    ! end (ln_bfrimp) 
    428  
    429                      
    430       !                                   !* d/dt(Ua), Ub, Vn (Vertical mean velocity) 
    431       !                                   ! --------------------------  
    432       zua(:,:) = zua(:,:) * hur(:,:) 
    433       zva(:,:) = zva(:,:) * hvr(:,:) 
    434       ! 
    435       ub_b(:,:) = ub_b(:,:) * zhu_b(:,:) 
    436       vb_b(:,:) = vb_b(:,:) * zhv_b(:,:) 
    437  
     439      END DO  
     440      ! 
     441      !                 ! Add bottom stress contribution from baroclinic velocities:       
     442      IF (ln_bt_fw) THEN 
     443         DO jj = 2, jpjm1                           
     444            DO ji = fs_2, fs_jpim1   ! vector opt. 
     445               ikbu = mbku(ji,jj)        
     446               ikbv = mbkv(ji,jj)     
     447               zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities 
     448               zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj) 
     449            END DO 
     450         END DO 
     451      ELSE 
     452         DO jj = 2, jpjm1 
     453            DO ji = fs_2, fs_jpim1   ! vector opt. 
     454               ikbu = mbku(ji,jj)        
     455               ikbv = mbkv(ji,jj)     
     456               zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities 
     457               zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj) 
     458            END DO 
     459         END DO 
     460      ENDIF 
     461      ! 
     462      ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 
     463      zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:) 
     464      zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:) 
     465      !        
     466      IF (ln_bt_fw) THEN                        ! Add wind forcing 
     467         zu_frc(:,:) =  zu_frc(:,:) + zraur * utau(:,:) * hur(:,:) 
     468         zv_frc(:,:) =  zv_frc(:,:) + zraur * vtau(:,:) * hvr(:,:) 
     469      ELSE 
     470         zu_frc(:,:) =  zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * hur(:,:) 
     471         zv_frc(:,:) =  zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:) 
     472      ENDIF   
     473      ! 
     474      IF ( ln_apr_dyn ) THEN                    ! Add atm pressure forcing 
     475         IF (ln_bt_fw) THEN 
     476            DO jj = 2, jpjm1               
     477               DO ji = fs_2, fs_jpim1   ! vector opt. 
     478                  zu_spg =  grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) /e1u(ji,jj) 
     479                  zv_spg =  grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) /e2v(ji,jj) 
     480                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 
     481                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg 
     482               END DO 
     483            END DO 
     484         ELSE 
     485            DO jj = 2, jpjm1               
     486               DO ji = fs_2, fs_jpim1   ! vector opt. 
     487                  zu_spg =  grav * z1_2 * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)    & 
     488                      &                    + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) /e1u(ji,jj) 
     489                  zv_spg =  grav * z1_2 * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj)    & 
     490                      &                    + ssh_ibb(ji  ,jj+1) - ssh_ibb(ji,jj)  ) /e2v(ji,jj) 
     491                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 
     492                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg 
     493               END DO 
     494            END DO 
     495         ENDIF  
     496      ENDIF 
     497      !                                   !* Right-Hand-Side of the barotropic ssh equation 
     498      !                                   ! ----------------------------------------------- 
     499      !                                         ! Surface net water flux and rivers 
     500      IF (ln_bt_fw) THEN 
     501         zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) ) 
     502      ELSE 
     503         zssh_frc(:,:) = zraur * z1_2 * (emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)) 
     504      ENDIF 
     505#if defined key_asminc 
     506      !                                         ! Include the IAU weighted SSH increment 
     507      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 
     508         zssh_frc(:,:) = zssh_frc(:,:) + ssh_iau(:,:) 
     509      ENDIF 
     510#endif 
     511      ! 
    438512      ! ----------------------------------------------------------------------- 
    439       !  Phase 2 : Integration of the barotropic equations with time splitting 
     513      !  Phase 2 : Integration of the barotropic equations  
    440514      ! ----------------------------------------------------------------------- 
    441515      ! 
     
    456530      ! Initialize depths: 
    457531      IF ( lk_vvl.AND.(.NOT.ln_bt_fw) ) THEN 
    458          hur_e  (:,:) = zhu_b  (:,:)     
    459          hvr_e  (:,:) = zhv_b  (:,:) 
    460          hu_e(:,:) = umask(:,:,1) / ( zhu_b(:,:) + 1. - umask(:,:,1) ) 
    461          hv_e(:,:) = vmask(:,:,1) / ( zhv_b(:,:) + 1. - vmask(:,:,1) ) 
    462       ELSE 
    463          hu_e   (:,:) = hu   (:,:)        
    464          hv_e   (:,:) = hv   (:,:)  
    465          hur_e  (:,:) = hur  (:,:)     
    466          hvr_e  (:,:) = hvr  (:,:) 
     532         hu_e  (:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1._wp - umask(:,:,1) ) 
     533         hv_e  (:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1._wp - vmask(:,:,1) ) 
     534         hur_e (:,:) = zhur_b(:,:) 
     535         hvr_e (:,:) = zhvr_b(:,:) 
     536      ELSE 
     537         hu_e  (:,:) = hu   (:,:)        
     538         hv_e  (:,:) = hv   (:,:)  
     539         hur_e (:,:) = hur  (:,:)     
     540         hvr_e (:,:) = hvr  (:,:) 
    467541      ENDIF 
    468542      ! 
     
    479553      zv_sum(:,:) = 0._wp 
    480554      !                                             ! ==================== ! 
    481        
    482       zsshb_e(:,:) = sshn_b(:,:)           ! sea surface height (before and now) 
    483       zub_e  (:,:) = un_b (:,:) 
    484       zvb_e  (:,:) = vn_b (:,:) 
    485  
    486       zu_sum  (:,:) = un_b (:,:)           ! summation 
    487       zv_sum  (:,:) = vn_b (:,:) 
    488       zssh_sum(:,:) = sshn (:,:) 
    489       ua_b  (:,:) = 0._wp                  ! After barotropic velocities (or transport if flux form)           
    490       va_b  (:,:) = 0._wp 
    491  
    492 #if defined key_obc 
    493       ! set ssh corrections to 0 
    494       ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop 
    495       IF( lp_obc_east  )   sshfoe_b(:,:) = 0._wp 
    496       IF( lp_obc_west  )   sshfow_b(:,:) = 0._wp 
    497       IF( lp_obc_south )   sshfos_b(:,:) = 0._wp 
    498       IF( lp_obc_north )   sshfon_b(:,:) = 0._wp 
    499 #endif 
    500  
    501       !                                             ! ==================== ! 
    502       DO jn = 1, icycle                             !  sub-time-step loop  ! (from NOW to AFTER+1) 
     555      DO jn = 1, icycle                             !  sub-time-step loop  ! 
    503556         !                                          ! ==================== ! 
    504          z2dt_e = 2. * ( rdt / nn_baro ) 
    505          IF( jn == 1 )   z2dt_e = rdt / nn_baro 
    506  
    507557         !                                                !* Update the forcing (BDY and tides) 
    508558         !                                                !  ------------------ 
    509          IF( lk_obc                    )  CALL obc_dta_bt( kt, jn   ) 
    510          IF ( lk_bdy .AND. lk_tide )      CALL bdy_dta_tides( kt, kit=jn, time_offset=(ioffset+1) ) 
    511          IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, koffset=ioffset ) 
     559         ! Update only tidal forcing at open boundaries 
     560#if defined key_tide 
     561         IF ( lk_bdy .AND. lk_tide )      CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) ) 
     562         IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, koffset=noffset ) 
     563#endif 
    512564         ! 
    513565         ! Set extrapolation coefficients for predictor step: 
     
    523575 
    524576         ! Extrapolate barotropic velocities at step jit+0.5: 
    525          ua_e(:,:) = za1 * zun_e(:,:) + za2 * zub_e(:,:) + za3 * ubb_e(:,:) 
    526          va_e(:,:) = za1 * zvn_e(:,:) + za2 * zvb_e(:,:) + za3 * vbb_e(:,:) 
     577         ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:) 
     578         va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) 
    527579 
    528580         IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only) 
    529581            !                                             !  ------------------ 
    530582            ! Extrapolate Sea Level at step jit+0.5: 
    531             zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * zsshb_e(:,:) + za3 * sshbb_e(:,:) 
     583            zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 
    532584            ! 
    533585            DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points 
     
    551603         ! considering fluxes below: 
    552604         ! 
    553          zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)   ! fluxes at jn+0.5  ! ACC INSTANT FAILURE IF THESE ARE USED (VERY LARGE NEGATIVE SALINITY AFTER 1 TS) 
     605         zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)         ! fluxes at jn+0.5 
    554606         zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 
    555 !        zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:)     ! fluxes at jn+0.5  ! ACC WORKS BETTER BUT FAILS AFTER O(100) TIMESTEPS 
    556 !        zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:) 
    557607         DO jj = 2, jpjm1                                  
    558608            DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    563613         END DO 
    564614         ! 
    565 #if defined key_obc 
    566          !                                                     ! OBC : zhdiv must be zero behind the open boundary 
    567 !!  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 
    568          IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1  ) = 0._wp      ! east 
    569          IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1  ) = 0._wp      ! west 
    570          IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0._wp      ! north 
    571          IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1  ) = 0._wp      ! south 
    572 #endif 
    573 #if defined key_bdy 
    574          zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:)               ! BDY mask 
    575 #endif 
    576          ! 
    577          DO jj = 2, jpjm1                                      ! leap-frog on ssh_e 
    578             DO ji = fs_2, fs_jpim1   ! vector opt. 
    579                ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * ( emp(ji,jj)-rnf(ji,jj) ) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1)  
    580             END DO 
    581          END DO 
     615         ! Sum over sub-time-steps to compute advective velocities 
     616         za2 = wgtbtp2(jn) 
     617         zu_sum  (:,:) = zu_sum  (:,:) + za2 * ua_e  (:,:) * zhup2_e  (:,:) 
     618         zv_sum  (:,:) = zv_sum  (:,:) + za2 * va_e  (:,:) * zhvp2_e  (:,:) 
     619         ! 
     620         ! Set next sea level: 
     621         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1) 
     622         CALL lbc_lnk( ssha_e, 'T',  1._wp ) 
     623 
    582624#if defined key_bdy 
    583625         ! Duplicate sea level across open boundaries (this is only cosmetic if lk_vvl=.false.) 
    584626         IF (lk_bdy) CALL bdy_ssh( ssha_e ) 
    585627#endif 
     628         !   
     629         ! Sea Surface Height at u-,v-points (vvl case only) 
     630         IF ( lk_vvl ) THEN                                 
     631            DO jj = 2, jpjm1 
     632               DO ji = 2, jpim1      ! NO Vector Opt. 
     633                  zsshu_a(ji,jj) = z1_2  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                 & 
     634                     &                                   * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha_e(ji  ,jj) & 
     635                     &                                     + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha_e(ji+1,jj) ) 
     636                  zsshv_a(ji,jj) = z1_2  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                 & 
     637                     &                                   * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha_e(ji,jj  ) & 
     638                     &                                     + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha_e(ji,jj+1) ) 
     639               END DO 
     640            END DO 
     641            CALL lbc_lnk( zsshu_a, 'U', 1._wp )   ;   CALL lbc_lnk( zsshv_a, 'V', 1._wp ) 
     642         ENDIF    
    586643         !                                  
    587644         ! Half-step back interpolation of SSH for surface pressure computation: 
     
    605662 
    606663         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:) & 
    607           &            + za2 *  zsshb_e(:,:) + za3 *  sshbb_e(:,:) 
     664          &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:) 
    608665 
    609666         ! 
     
    635692            DO jj = 2, jpjm1 
    636693               DO ji = fs_2, fs_jpim1   ! vector opt. 
    637                   ! surface pressure gradient 
    638                   IF( lk_vvl) THEN 
    639                      zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    & 
    640                         &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj) 
    641                      zv_spg = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    & 
    642                         &              - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj) 
    643                   ELSE 
    644                      zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 
    645                      zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
    646                   ENDIF 
    647                   ! add tidal astronomical forcing 
    648                   IF ( ln_tide_pot .AND. lk_tide ) THEN  
    649                   zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 
    650                   zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 
    651                   ENDIF 
    652                   ! energy conserving formulation for planetary vorticity term 
    653694                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj) 
    654695                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
    655696                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj) 
    656697                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
    657                   zu_cor = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj) 
    658                   zv_cor =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj) 
    659                   ! after velocities with implicit bottom friction 
    660  
    661                   IF( ln_bfrimp ) THEN      ! implicit bottom friction 
    662                      !   A new method to implement the implicit bottom friction.  
    663                      !   H. Liu 
    664                      !   Sept 2011 
    665                      ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) +                                            & 
    666                       &                               z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )            & 
    667                       &                               / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 
    668                      ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)    
    669                      ! 
    670                      va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) +                                            & 
    671                       &                               z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )            & 
    672                       &                               / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 
    673                      va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)    
    674                      ! 
    675                   ELSE 
    676                      ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1)   & 
    677                       &           / ( 1._wp         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
    678                      va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1)   & 
    679                       &           / ( 1._wp         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
    680                   ENDIF 
     698                  zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
     699                  zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
    681700               END DO 
    682701            END DO 
     
    685704            DO jj = 2, jpjm1 
    686705               DO ji = fs_2, fs_jpim1   ! vector opt. 
    687                    ! surface pressure gradient 
    688                   IF( lk_vvl) THEN 
    689                      zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    & 
    690                         &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj) 
    691                      zv_spg = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    & 
    692                         &              - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj) 
    693                   ELSE 
    694                      zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 
    695                      zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
    696                   ENDIF 
    697                   ! add tidal astronomical forcing 
    698                   IF ( ln_tide_pot .AND. lk_tide ) THEN 
    699                   zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 
    700                   zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 
    701                   ENDIF 
    702                   ! enstrophy conserving formulation for planetary vorticity term 
    703                   zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
    704                   zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
    705                   zu_cor  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj) 
    706                   zv_cor  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) * hvr_e(ji,jj) 
    707                   ! after velocities with implicit bottom friction 
    708                   IF( ln_bfrimp ) THEN 
    709                      !   A new method to implement the implicit bottom friction.  
    710                      !   H. Liu 
    711                      !   Sept 2011 
    712                      ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) +                                            & 
    713                       &                               z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )            & 
    714                       &                               / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 
    715                      ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)    
    716                      ! 
    717                      va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) +                                            & 
    718                       &                               z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )            & 
    719                       &                               / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 
    720                      va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)    
    721                      ! 
    722                   ELSE 
    723                      ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1)   & 
    724                      &            / ( 1._wp        - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
    725                      va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1)   & 
    726                      &            / ( 1._wp        - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
    727                   ENDIF 
     706                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) & 
     707                   &             + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
     708                  zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) & 
     709                   &             + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
     710                  zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) ) 
     711                  zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) ) 
    728712               END DO 
    729713            END DO 
     
    732716            DO jj = 2, jpjm1 
    733717               DO ji = fs_2, fs_jpim1   ! vector opt. 
    734                   ! surface pressure gradient 
    735                   IF( lk_vvl) THEN 
    736                      zu_spg = -grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn_e(ji+1,jj  )    & 
    737                         &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj) 
    738                      zv_spg = -grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    & 
    739                         &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj) 
    740                   ELSE 
    741                      zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 
    742                      zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
    743                   ENDIF 
    744                   ! add tidal astronomical forcing 
    745                   IF ( ln_tide_pot .AND. lk_tide ) THEN 
    746                   zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 
    747                   zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 
    748                   ENDIF 
    749                   ! energy/enstrophy conserving formulation for planetary vorticity term 
    750                   zu_cor = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   & 
    751                      &                           + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) * hur_e(ji,jj) 
    752                   zv_cor = - z1_4 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   & 
    753                      &                           + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) * hvr_e(ji,jj) 
    754                   ! after velocities with implicit bottom friction 
    755                   IF( ln_bfrimp ) THEN 
    756                      !   A new method to implement the implicit bottom friction.  
    757                      !   H. Liu 
    758                      !   Sept 2011 
    759                      ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) +                                            & 
    760                       &                               z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )            & 
    761                       &                               / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 
    762                      ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)    
    763                      ! 
    764                      va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) +                                            & 
    765                       &                               z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )            & 
    766                       &                               / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 
    767                      va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)    
    768                      ! 
    769                   ELSE 
    770                      ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1)   & 
    771                      &            / ( 1._wp        - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
    772                      va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1)   & 
    773                      &            / ( 1._wp        - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
    774                   ENDIF 
     718                  zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) & 
     719                     &                                    + ftnw(ji+1,jj) * zwy(ji+1,jj  ) & 
     720                     &                                    + ftse(ji,jj  ) * zwy(ji  ,jj-1) &  
     721                     &                                    + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 
     722                  zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) &  
     723                     &                                    + ftse(ji,jj+1) * zwx(ji  ,jj+1) & 
     724                     &                                    + ftnw(ji,jj  ) * zwx(ji-1,jj  ) &  
     725                     &                                    + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) 
    775726               END DO 
    776727            END DO 
    777728            !  
    778729         ENDIF 
    779          !                                                     !* domain lateral boundary 
    780          !                                                     !  ----------------------- 
    781  
    782                                                                ! OBC open boundaries 
    783          IF( lk_obc               )   CALL obc_fla_ts ( ua_e, va_e, sshn_e, ssha_e ) 
    784  
    785                                                                ! BDY open boundaries 
    786 #if defined key_bdy 
    787          pssh => sshn_e 
     730         ! 
     731         ! Add tidal astronomical forcing if defined 
     732         IF ( lk_tide.AND.ln_tide_pot ) THEN 
     733            DO jj = 2, jpjm1 
     734               DO ji = fs_2, fs_jpim1   ! vector opt. 
     735                  zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 
     736                  zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 
     737                  zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg 
     738                  zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg 
     739               END DO 
     740            END DO 
     741         ENDIF 
     742         ! 
     743         ! Add bottom stresses: 
     744         zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:) 
     745         zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:) 
     746         ! 
     747         ! Surface pressure trend: 
     748         DO jj = 2, jpjm1 
     749            DO ji = fs_2, fs_jpim1   ! vector opt. 
     750               ! Add surface pressure gradient 
     751               zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 
     752               zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 
     753               zwx(ji,jj) = zu_spg 
     754               zwy(ji,jj) = zv_spg 
     755            END DO 
     756         END DO 
     757         ! 
     758         ! Set next velocities: 
     759         IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN    ! Vector form 
     760            DO jj = 2, jpjm1 
     761               DO ji = fs_2, fs_jpim1   ! vector opt. 
     762                  ua_e(ji,jj) = (                                zun_e(ji,jj)   &  
     763                            &     + rdtbt * (                      zwx(ji,jj)   & 
     764                            &                                 + zu_trd(ji,jj)   & 
     765                            &                                 + zu_frc(ji,jj) ) &  
     766                            &   ) * umask(ji,jj,1) 
     767 
     768                  va_e(ji,jj) = (                                zvn_e(ji,jj)   & 
     769                            &     + rdtbt * (                      zwy(ji,jj)   & 
     770                            &                                 + zv_trd(ji,jj)   & 
     771                            &                                 + zv_frc(ji,jj) ) & 
     772                            &   ) * vmask(ji,jj,1) 
     773               END DO 
     774            END DO 
     775 
     776         ELSE                 ! Flux form 
     777            DO jj = 2, jpjm1 
     778               DO ji = fs_2, fs_jpim1   ! vector opt. 
     779 
     780                  zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1)) 
     781                  zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 
     782 
     783                  ua_e(ji,jj) = (                hu_e(ji,jj)  *  zun_e(ji,jj)   &  
     784                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   &  
     785                            &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   & 
     786                            &               +      hu(ji,jj)  * zu_frc(ji,jj) ) & 
     787                            &   ) * zhura 
     788 
     789                  va_e(ji,jj) = (                hv_e(ji,jj)  *  zvn_e(ji,jj)   & 
     790                            &     + rdtbt * ( zhvst_e(ji,jj)  *    zwy(ji,jj)   & 
     791                            &               + zhvp2_e(ji,jj)  * zv_trd(ji,jj)   & 
     792                            &               +      hv(ji,jj)  * zv_frc(ji,jj) ) & 
     793                            &   ) * zhvra 
     794               END DO 
     795            END DO 
     796         ENDIF 
     797         ! 
     798         IF( lk_vvl ) THEN                             !* Update ocean depth (variable volume case only) 
     799            !                                          !  ----------------------------------------------         
     800            hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 
     801            hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 
     802            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 
     803            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) 
     804            ! 
     805         ENDIF 
     806         !                                                 !* domain lateral boundary 
     807         !                                                 !  ----------------------- 
     808         ! 
     809         CALL lbc_lnk( ua_e  , 'U', -1._wp )               ! local domain boundaries  
     810         CALL lbc_lnk( va_e  , 'V', -1._wp ) 
     811 
     812#if defined key_bdy   
     813 
     814         pssh => ssha_e 
    788815         phur => hur_e 
    789816         phvr => hvr_e 
    790817         pu2d => ua_e 
    791818         pv2d => va_e 
    792  
    793          IF( lk_bdy )   CALL bdy_dyn2d( kt )  
     819                                       
     820         IF( lk_bdy )   CALL bdy_dyn2d( kt )               ! open boundaries 
    794821#endif 
    795  
    796          ! 
    797          CALL lbc_lnk( ua_e  , 'U', -1. )                      ! local domain boundaries  
    798          CALL lbc_lnk( va_e  , 'V', -1. ) 
    799          CALL lbc_lnk( ssha_e, 'T',  1. ) 
    800  
    801          zu_sum  (:,:) = zu_sum  (:,:) + ua_e  (:,:)           ! Sum over sub-time-steps 
    802          zv_sum  (:,:) = zv_sum  (:,:) + va_e  (:,:)  
    803          zssh_sum(:,:) = zssh_sum(:,:) + ssha_e(:,:)  
    804  
    805          !                                                !* Time filter and swap 
    806          !                                                !  -------------------- 
    807          sshbb_e(:,:) = zsshb_e(:,:) 
    808          ubb_e  (:,:) = zub_e  (:,:) 
    809          vbb_e  (:,:) = zvb_e  (:,:) 
    810          IF( jn == 1 ) THEN                                     ! Swap only (1st Euler time step) 
    811             zsshb_e(:,:) = sshn_e(:,:) 
    812             zub_e  (:,:) = zun_e (:,:) 
    813             zvb_e  (:,:) = zvn_e (:,:) 
    814             sshn_e (:,:) = ssha_e(:,:) 
    815             zun_e  (:,:) = ua_e  (:,:) 
    816             zvn_e  (:,:) = va_e  (:,:) 
    817          ELSE                                                   ! Swap + Filter 
    818             zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:) 
    819             zub_e  (:,:) = atfp * ( zub_e  (:,:) + ua_e  (:,:) ) + atfp1 * zun_e (:,:) 
    820             zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + va_e  (:,:) ) + atfp1 * zvn_e (:,:) 
    821             sshn_e (:,:) = ssha_e(:,:) 
    822             zun_e  (:,:) = ua_e  (:,:) 
    823             zvn_e  (:,:) = va_e  (:,:) 
    824          ENDIF 
    825  
    826          za1 = wgtbtp1(jn) 
     822         !                                             !* Swap 
     823         !                                             !  ---- 
     824         ubb_e  (:,:) = ub_e  (:,:) 
     825         ub_e   (:,:) = zun_e (:,:) 
     826         zun_e  (:,:) = ua_e  (:,:) 
     827         ! 
     828         vbb_e  (:,:) = vb_e  (:,:) 
     829         vb_e   (:,:) = zvn_e (:,:) 
     830         zvn_e  (:,:) = va_e  (:,:) 
     831         ! 
     832         sshbb_e(:,:) = sshb_e(:,:) 
     833         sshb_e (:,:) = sshn_e(:,:) 
     834         sshn_e (:,:) = ssha_e(:,:) 
     835 
     836         !                                             !* Sum over whole bt loop 
     837         !                                             !  ---------------------- 
     838         za1 = wgtbtp1(jn)                                     
    827839         IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN    ! Sum velocities 
    828             ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) 
    829             va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) 
    830          ELSE                                              ! Sum transports 
    831             ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) 
    832             va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) 
    833 !           ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:) 
    834 !           va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:) 
    835          ENDIF 
     840            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:)  
     841            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:)  
     842         ELSE                                ! Sum transports 
     843            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:) 
     844            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:) 
     845         ENDIF 
     846         !                                   ! Sum sea level 
    836847         ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:) 
    837  
    838          IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only) 
    839             !                                             !  ------------------ 
    840             DO jj = 1, jpjm1                                    ! Sea Surface Height at u- & v-points 
    841                DO ji = 1, fs_jpim1   ! Vector opt. 
    842                   zsshun_e(ji,jj) = 0.5_wp * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       & 
    843                      &                     * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn_e(ji  ,jj)    & 
    844                      &                     +   e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) ) 
    845                   zsshvn_e(ji,jj) = 0.5_wp * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       & 
    846                      &                     * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn_e(ji,jj  )    & 
    847                      &                     +   e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) ) 
    848                END DO 
    849             END DO 
    850             CALL lbc_lnk( zsshun_e, 'U', 1. )                   ! lateral boundaries conditions 
    851             CALL lbc_lnk( zsshvn_e, 'V', 1. )  
    852             ! 
    853             hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:)              ! Ocean depth at U- and V-points 
    854             hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:) 
    855             hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 
    856             hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) 
    857             ! 
    858          ENDIF 
    859848         !                                                 ! ==================== ! 
    860849      END DO                                               !        end loop      ! 
    861850      !                                                    ! ==================== ! 
    862  
    863 #if defined key_obc 
    864       IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:)     !!gm totally useless ????? 
    865       IF( lp_obc_west  )   sshfow_b(:,:) = zcoef * sshfow_b(:,:) 
    866       IF( lp_obc_north )   sshfon_b(:,:) = zcoef * sshfon_b(:,:) 
    867       IF( lp_obc_south )   sshfos_b(:,:) = zcoef * sshfos_b(:,:) 
    868 #endif 
    869  
    870851      ! ----------------------------------------------------------------------------- 
    871852      ! Phase 3. update the general trend with the barotropic trend 
    872853      ! ----------------------------------------------------------------------------- 
    873854      ! 
    874       !                                   !* Time average ==> after barotropic u, v, ssh 
    875       zcoef =  1._wp / ( icycle + 1 )  
    876       zu_sum(:,:) = zcoef * zu_sum  (:,:)  
    877       zv_sum(:,:) = zcoef * zv_sum  (:,:)  
    878       !  
     855      ! At this stage ssha holds a time averaged value 
     856      !                                                ! Sea Surface Height at u-,v- and f-points 
     857      IF( lk_vvl ) THEN                                ! (required only in key_vvl case) 
     858         DO jj = 1, jpjm1 
     859            DO ji = 1, jpim1      ! NO Vector Opt. 
     860               zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   & 
     861                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     & 
     862                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) ) 
     863               zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   & 
     864                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     & 
     865                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) 
     866            END DO 
     867         END DO 
     868         CALL lbc_lnk( zsshu_a, 'U', 1._wp )   ;   CALL lbc_lnk( zsshv_a, 'V', 1._wp ) ! Boundary conditions 
     869      ENDIF 
     870      ! 
    879871      ! Set advection velocity correction: 
    880       IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN 
    881          un_adv(:,:) = zu_sum(:,:) 
    882          vn_adv(:,:) = zv_sum(:,:) 
    883       ELSE 
    884          un_adv(:,:) = 0.5_wp * ( ub_b(:,:) + zu_sum(:,:)) 
    885          vn_adv(:,:) = 0.5_wp * ( vb_b(:,:) + zv_sum(:,:)) 
     872      IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN      
     873         un_adv(:,:) = zu_sum(:,:)*hur(:,:) 
     874         vn_adv(:,:) = zv_sum(:,:)*hvr(:,:) 
     875      ELSE 
     876         un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:)) * hur(:,:) 
     877         vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:)) * hvr(:,:) 
    886878      END IF 
    887       !                                   !* update the general momentum trend 
    888       DO jk=1,jpkm1 
    889          ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) * z1_2dt_b 
    890          va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) * z1_2dt_b 
    891       END DO 
    892       un_b  (:,:) =  zu_sum(:,:)  
    893       vn_b  (:,:) =  zv_sum(:,:)  
    894       sshn_b(:,:) = zcoef * zssh_sum(:,:)  
     879 
     880      IF (ln_bt_fw) THEN ! Save integrated transport for next computation 
     881         ub2_b(:,:) = zu_sum(:,:) 
     882         vb2_b(:,:) = zv_sum(:,:) 
     883      ENDIF 
     884      ! 
     885      ! Update barotropic trend: 
     886      IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN 
     887         DO jk=1,jpkm1 
     888            ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 
     889            va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 
     890         END DO 
     891      ELSE 
     892         hu_e  (:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1._wp - umask(:,:,1) ) 
     893         hv_e  (:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1._wp - vmask(:,:,1) ) 
     894         DO jk=1,jpkm1 
     895            ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_e(:,:) ) * z1_2dt_b 
     896            va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_e(:,:) ) * z1_2dt_b 
     897         END DO 
     898         ! Save barotropic velocities not transport: 
     899         ua_b  (:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) ) 
     900         va_b  (:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) ) 
     901      ENDIF 
    895902      ! 
    896903      !                                   !* write time-spliting arrays in the restart 
    897       IF( lrst_oce )   CALL ts_rst( kt, 'WRITE' ) 
    898       ! 
    899       CALL wrk_dealloc( jpi, jpj, zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv     ) 
    900       CALL wrk_dealloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e   ) 
    901       CALL wrk_dealloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum ) 
    902       CALL wrk_dealloc( jpi, jpj, zhu_b, zhv_b                                     ) 
    903       CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zsshp2_e                       ) 
    904       CALL wrk_dealloc( jpi, jpj, zhust_e, zhvst_e                                 ) 
     904      IF(lrst_oce .AND.ln_bt_fw)   CALL ts_rst( kt, 'WRITE' ) 
     905      ! 
     906      CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv ) 
     907      CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e ) 
     908      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc ) 
     909      CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 
     910      CALL wrk_dealloc( jpi, jpj, zhur_b, zhvr_b                                     ) 
     911      CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a                                   ) 
     912      CALL wrk_dealloc( jpi, jpj, zht, zhf ) 
    905913      ! 
    906914      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts') 
     
    9941002      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag 
    9951003      ! 
    996       REAL(wp), POINTER, DIMENSION(:,:) :: zzhu_b, zzhv_b 
    997       INTEGER ::  ji, jk        ! dummy loop indices 
    9981004      !!---------------------------------------------------------------------- 
    9991005      ! 
    10001006      IF( TRIM(cdrw) == 'READ' ) THEN 
    1001          IF( iom_varid( numror, 'un_b', ldstop = .FALSE. ) > 0 ) THEN 
    1002             CALL iom_get( numror, jpdom_autoglo, 'un_b'  , un_b  (:,:) )   ! external velocity issued 
    1003             CALL iom_get( numror, jpdom_autoglo, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop 
    1004          ELSE 
    1005             un_b (:,:) = 0._wp 
    1006             vn_b (:,:) = 0._wp 
    1007             ! vertical sum 
    1008             IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll 
    1009                DO jk = 1, jpkm1 
    1010                   DO ji = 1, jpij 
    1011                      un_b(ji,1) = un_b(ji,1) + fse3u_n(ji,1,jk) * un(ji,1,jk) 
    1012                      vn_b(ji,1) = vn_b(ji,1) + fse3v_n(ji,1,jk) * vn(ji,1,jk) 
    1013                   END DO 
    1014                END DO 
    1015             ELSE                             ! No  vector opt. 
    1016                DO jk = 1, jpkm1 
    1017                   un_b(:,:) = un_b(:,:) + fse3u_n(:,:,jk) * un(:,:,jk) 
    1018                   vn_b(:,:) = vn_b(:,:) + fse3v_n(:,:,jk) * vn(:,:,jk) 
    1019                END DO 
    1020             ENDIF 
    1021             un_b (:,:) = un_b(:,:) * hur(:,:) 
    1022             vn_b (:,:) = vn_b(:,:) * hvr(:,:) 
    1023          ENDIF 
    1024  
    1025          ! Vertically integrated velocity (before) 
    1026          IF (neuler/=0) THEN 
    1027             ub_b (:,:) = 0._wp 
    1028             vb_b (:,:) = 0._wp 
    1029  
    1030             ! vertical sum 
    1031             IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll 
    1032                DO jk = 1, jpkm1 
    1033                   DO ji = 1, jpij 
    1034                      ub_b(ji,1) = ub_b(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk) 
    1035                      vb_b(ji,1) = vb_b(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk) 
    1036                   END DO 
    1037                END DO 
    1038             ELSE                             ! No  vector opt. 
    1039                DO jk = 1, jpkm1 
    1040                   ub_b(:,:) = ub_b(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) 
    1041                   vb_b(:,:) = vb_b(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) 
    1042                END DO 
    1043             ENDIF 
    1044  
    1045             IF( lk_vvl ) THEN 
    1046                CALL wrk_alloc( jpi, jpj, zzhu_b, zzhv_b ) 
    1047                ub_b (:,:) = 0. 
    1048                vb_b (:,:) = 0. 
    1049                zzhu_b(:,:) = 0. 
    1050                zzhv_b(:,:) = 0. 
    1051                ! vertical sum 
    1052                IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll 
    1053                   DO jk = 1, jpkm1 
    1054                      DO ji = 1, jpij 
    1055                         ub_b  (ji,1) = ub_b (ji,1)  + fse3u_b(ji,1,jk) * ub   (ji,1,jk) 
    1056                         vb_b  (ji,1) = vb_b (ji,1)  + fse3v_b(ji,1,jk) * vb   (ji,1,jk) 
    1057                         zzhu_b(ji,1) = zzhu_b(ji,1) + fse3u_b(ji,1,jk) * umask(ji,1,jk) 
    1058                         zzhv_b(ji,1) = zzhv_b(ji,1) + fse3v_b(ji,1,jk) * vmask(ji,1,jk) 
    1059                      END DO 
    1060                   END DO 
    1061                ELSE                             ! No  vector opt. 
    1062                   DO jk = 1, jpkm1 
    1063                      ub_b  (:,:) = ub_b  (:,:) + fse3u_b(:,:,jk) * ub   (:,:,jk) 
    1064                      vb_b  (:,:) = vb_b  (:,:) + fse3v_b(:,:,jk) * vb   (:,:,jk) 
    1065                      zzhu_b(:,:) = zzhu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk) 
    1066                      zzhv_b(:,:) = zzhv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 
    1067                   END DO 
    1068                ENDIF 
    1069                ub_b(:,:) = ub_b(:,:) / ( zzhu_b(:,:) + 1. - umask(:,:,1) ) 
    1070                vb_b(:,:) = vb_b(:,:) / ( zzhv_b(:,:) + 1. - vmask(:,:,1) ) 
    1071                CALL wrk_dealloc( jpi, jpj, zzhu_b, zzhv_b ) 
    1072             ELSE              
    1073                ub_b (:,:) = 0.e0 
    1074                vb_b (:,:) = 0.e0 
    1075                ! vertical sum 
    1076                IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll 
    1077                   DO jk = 1, jpkm1 
    1078                      DO ji = 1, jpij 
    1079                         ub_b(ji,1) = ub_b(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk) 
    1080                         vb_b(ji,1) = vb_b(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk) 
    1081                      END DO 
    1082                   END DO 
    1083                ELSE                             ! No  vector opt. 
    1084                   DO jk = 1, jpkm1 
    1085                      ub_b(:,:) = ub_b(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) 
    1086                      vb_b(:,:) = vb_b(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) 
    1087                   END DO 
    1088                ENDIF 
    1089                ub_b(:,:) = ub_b(:,:) * hur(:,:) 
    1090                vb_b(:,:) = vb_b(:,:) * hvr(:,:) 
    1091             ENDIF 
    1092          ELSE                                 ! neuler==0 
    1093             ub_b (:,:) = un_b (:,:) 
    1094             vb_b (:,:) = vn_b (:,:) 
    1095          ENDIF 
    1096  
    1097          IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN 
    1098             CALL iom_get( numror, jpdom_autoglo, 'sshn_b' , sshn_b (:,:) )   ! filtered ssh 
    1099          ELSE 
    1100             sshn_b(:,:) = sshb(:,:)   ! if not in restart set previous time mean to current baroclinic before value    
    1101          ENDIF  
    1102  
     1007         CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:) )    
     1008         CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:) )  
    11031009         IF( .NOT.ln_bt_av .AND. iom_varid( numror, 'sshbb_e', ldstop = .FALSE. ) > 0) THEN 
    11041010            CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:) )    
     
    11161022            vb_e    = vb_b 
    11171023         ENDIF 
     1024      ! 
    11181025      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 
    1119          CALL iom_rstput( kt, nitrst, numrow, 'un_b'   , un_b  (:,:) )   ! external velocity and ssh 
    1120          CALL iom_rstput( kt, nitrst, numrow, 'vn_b'   , vn_b  (:,:) )   ! issued from barotropic loop 
    1121          CALL iom_rstput( kt, nitrst, numrow, 'sshn_b' , sshn_b(:,:) )   !  
     1026         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:) ) 
     1027         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:) ) 
     1028         ! 
    11221029         IF (.NOT.ln_bt_av) THEN 
    11231030            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:) )  
     
    11531060      CALL wrk_alloc( jpi, jpj, zcu ) 
    11541061      ! 
    1155       IF (lk_vvl) THEN 
    1156          DO jj = 1, jpj 
    1157             DO ji =1, jpi  
    1158                zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj)) 
    1159                zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj)) 
    1160                zcu(ji,jj) = sqrt(grav*ht_0(ji,jj)*(zxr2 + zyr2) ) 
    1161             END DO 
    1162          END DO 
    1163       ELSE 
    1164          DO jj = 1, jpj 
    1165             DO ji =1, jpi  
    1166                zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj)) 
    1167                zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj)) 
    1168                zcu(ji,jj) = sqrt( grav*ht_0(ji,jj)*(zxr2 + zyr2) ) 
    1169             END DO 
    1170          END DO 
    1171       ENDIF 
     1062      DO jj = 1, jpj 
     1063         DO ji =1, jpi  
     1064            zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj)) 
     1065            zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj)) 
     1066            zcu(ji,jj) = sqrt(grav*ht_0(ji,jj)*(zxr2 + zyr2) ) 
     1067         END DO 
     1068      END DO 
    11721069 
    11731070      zcmax = MAXVAL(zcu(:,:)) 
     
    12471144      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag 
    12481145      WRITE(*,*) 'ts_rst    : You should not have seen this print! error?', kt, cdrw 
    1249    END SUBROUTINE ts_rst     
     1146   END SUBROUTINE ts_rst   
    12501147   SUBROUTINE dyn_spg_ts_init( kt )       ! Empty routine 
    12511148      INTEGER         , INTENT(in) ::   kt         ! ocean time-step 
Note: See TracChangeset for help on using the changeset viewer.