Ignore:
Timestamp:
2019-09-19T11:18:03+02:00 (13 months ago)
Author:
jchanut
Message:

#2222, merged with trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DYN/dynspg_ts.F90

    r10742 r11573  
    6464   USE diatmb          ! Top,middle,bottom output 
    6565 
     66   USE iom   ! to remove 
     67 
    6668   IMPLICIT NONE 
    6769   PRIVATE 
     
    104106      ! 
    105107      ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) ) 
    106       ! 
    107108      IF( ln_dynvor_een .OR. ln_dynvor_eeT )   & 
    108          &     ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &  
    109          &               ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) ) 
     109         &     ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , ftsw(jpi,jpj) , ftse(jpi,jpj), STAT=ierr(2)   ) 
    110110         ! 
    111111      ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj)                    , STAT=ierr(3) ) 
     
    149149      LOGICAL  ::   ll_fw_start           ! =T : forward integration  
    150150      LOGICAL  ::   ll_init               ! =T : special startup of 2d equations 
    151       LOGICAL  ::   ll_tmp1, ll_tmp2      ! local logical variables used in W/D 
    152       INTEGER  ::   ikbu, iktu, noffset   ! local integers 
    153       INTEGER  ::   ikbv, iktv            !   -      - 
    154       REAL(wp) ::   r1_2dt_b, z2dt_bf               ! local scalars 
    155       REAL(wp) ::   zx1, zx2, zu_spg, zhura, z1_hu  !   -      - 
    156       REAL(wp) ::   zy1, zy2, zv_spg, zhvra, z1_hv  !   -      - 
     151      INTEGER  ::   noffset               ! local integers  : time offset for bdy update 
     152      REAL(wp) ::   r1_2dt_b, z1_hu, z1_hv          ! local scalars 
    157153      REAL(wp) ::   za0, za1, za2, za3              !   -      - 
    158       REAL(wp) ::   zmdi, zztmp            , z1_ht  !   -      - 
    159       REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e, zhf 
    160       REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zssh_frc 
    161       REAL(wp), DIMENSION(jpi,jpj) :: zwy, zv_trd, zv_frc, zhdiv 
    162       REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhust_e, zhtp2_e 
    163       REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zhvst_e 
     154      REAL(wp) ::   zmdi, zztmp, zldg               !   -      - 
     155      REAL(wp) ::   zhu_bck, zhv_bck, zhdiv         !   -      - 
     156      REAL(wp) ::   zun_save, zvn_save              !   -      - 
     157      REAL(wp), DIMENSION(jpi,jpj) :: zu_trd, zu_frc, zu_spg, zssh_frc 
     158      REAL(wp), DIMENSION(jpi,jpj) :: zv_trd, zv_frc, zv_spg 
     159      REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhtp2_e 
     160      REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zsshp2_e 
    164161      REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v   ! top/bottom stress at u- & v-points 
     162      REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV         ! fluxes 
    165163      ! 
    166164      REAL(wp) ::   zwdramp                     ! local scalar - only used if ln_wd_dl = .True.  
     
    182180      zwdramp = r_rn_wdmin1               ! simplest ramp  
    183181!     zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp 
    184       !                                         ! reciprocal of baroclinic time step  
    185       IF( kt == nit000 .AND. neuler == 0 ) THEN   ;   z2dt_bf =          rdt 
    186       ELSE                                        ;   z2dt_bf = 2.0_wp * rdt 
    187       ENDIF 
    188       r1_2dt_b = 1.0_wp / z2dt_bf  
     182      !                                         ! inverse of baroclinic time step  
     183      IF( kt == nit000 .AND. neuler == 0 ) THEN   ;   r1_2dt_b = 1._wp / (         rdt ) 
     184      ELSE                                        ;   r1_2dt_b = 1._wp / ( 2._wp * rdt ) 
     185      ENDIF 
    189186      ! 
    190187      ll_init     = ln_bt_av                    ! if no time averaging, then no specific restart  
     
    210207            ll_fw_start =.FALSE. 
    211208         ENDIF 
    212          ! 
    213          ! Set averaging weights and cycle length: 
     209         !                    ! Set averaging weights and cycle length: 
    214210         CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 
    215211         ! 
    216       ENDIF 
    217       ! 
    218       IF( ln_isfcav ) THEN    ! top+bottom friction (ocean cavities) 
    219          DO jj = 2, jpjm1 
    220             DO ji = fs_2, fs_jpim1   ! vector opt. 
    221                zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) 
    222                zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) 
    223             END DO 
    224          END DO 
    225       ELSE                    ! bottom friction only 
    226          DO jj = 2, jpjm1 
    227             DO ji = fs_2, fs_jpim1   ! vector opt. 
    228                zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) 
    229                zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) 
    230             END DO 
    231          END DO 
    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 the 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 zhf resp.) at each barotropic time step 
    240       ! 
    241       IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN 
    242          ! 
    243          SELECT CASE( nvor_scheme ) 
    244          CASE( np_EEN )                != EEN scheme using e3f (energy & enstrophy scheme) 
    245             SELECT CASE( nn_een_e3f )              !* ff_f/e3 at F-point 
    246             CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4) 
    247                DO jj = 1, jpjm1 
    248                   DO ji = 1, jpim1 
    249                      zwz(ji,jj) =   ( ht_n(ji  ,jj+1) + ht_n(ji+1,jj+1) +                    & 
    250                         &             ht_n(ji  ,jj  ) + ht_n(ji+1,jj  )   ) * 0.25_wp   
    251                      IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 
    252                   END DO 
    253                END DO 
    254             CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask) 
    255                DO jj = 1, jpjm1 
    256                   DO ji = 1, jpim1 
    257                      zwz(ji,jj) =             (  ht_n  (ji  ,jj+1) + ht_n  (ji+1,jj+1)      & 
    258                         &                      + ht_n  (ji  ,jj  ) + ht_n  (ji+1,jj  )  )   & 
    259                         &       / ( MAX( 1._wp,  ssmask(ji  ,jj+1) + ssmask(ji+1,jj+1)      & 
    260                         &                      + ssmask(ji  ,jj  ) + ssmask(ji+1,jj  )  )   ) 
    261                      IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 
    262                   END DO 
    263                END DO 
    264             END SELECT 
    265             CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp ) 
    266             ! 
    267             ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 
    268             DO jj = 2, jpj 
    269                DO ji = 2, jpi 
    270                   ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1) 
    271                   ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  ) 
    272                   ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1) 
    273                   ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  ) 
    274                END DO 
    275             END DO 
    276             ! 
    277          CASE( np_EET )                  != EEN scheme using e3t (energy conserving scheme) 
    278             ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 
    279             DO jj = 2, jpj 
    280                DO ji = 2, jpi 
    281                   z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) ) 
    282                   ftne(ji,jj) = ( ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) ) * z1_ht 
    283                   ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) ) * z1_ht 
    284                   ftse(ji,jj) = ( ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht 
    285                   ftsw(ji,jj) = ( ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) ) * z1_ht 
    286                END DO 
    287             END DO 
    288             ! 
    289          CASE( np_ENE, np_ENS , np_MIX )  != all other schemes (ENE, ENS, MIX) except ENT ! 
    290             ! 
    291             zwz(:,:) = 0._wp 
    292             zhf(:,:) = 0._wp 
    293              
    294 !!gm  assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed  
    295 !!gm    A priori a better value should be something like : 
    296 !!gm          zhf(i,j) = masked sum of  ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)  
    297 !!gm                     divided by the sum of the corresponding mask  
    298 !!gm  
    299 !!             
    300             IF( .NOT.ln_sco ) THEN 
    301    
    302    !!gm  agree the JC comment  : this should be done in a much clear way 
    303    
    304    ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 
    305    !     Set it to zero for the time being  
    306    !              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level 
    307    !              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth 
    308    !              ENDIF 
    309    !              zhf(:,:) = gdepw_0(:,:,jk+1) 
    310                ! 
    311             ELSE 
    312                ! 
    313                !zhf(:,:) = hbatf(:,:) 
    314                DO jj = 1, jpjm1 
    315                   DO ji = 1, jpim1 
    316                      zhf(ji,jj) =    (   ht_0  (ji,jj  ) + ht_0  (ji+1,jj  )          & 
    317                         &              + ht_0  (ji,jj+1) + ht_0  (ji+1,jj+1)   )      & 
    318                         &       / MAX(   ssmask(ji,jj  ) + ssmask(ji+1,jj  )          & 
    319                         &              + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp  ) 
    320                   END DO 
    321                END DO 
    322             ENDIF 
    323             ! 
    324             DO jj = 1, jpjm1 
    325                zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 
    326             END DO 
    327             ! 
    328             DO jk = 1, jpkm1 
    329                DO jj = 1, jpjm1 
    330                   zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 
    331                END DO 
    332             END DO 
    333             CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp ) 
    334             ! JC: TBC. hf should be greater than 0  
    335             DO jj = 1, jpj 
    336                DO ji = 1, jpi 
    337                   IF( zhf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zhf(ji,jj) ! zhf is actually hf here but it saves an array 
    338                END DO 
    339             END DO 
    340             zwz(:,:) = ff_f(:,:) * zwz(:,:) 
    341          END SELECT 
    342212      ENDIF 
    343213      ! 
     
    348218         CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 
    349219      ENDIF 
     220      ! 
    350221                           
    351222      ! ----------------------------------------------------------------------------- 
     
    354225      !       
    355226      ! 
    356       !                                   !* e3*d/dt(Ua) (Vertically integrated) 
    357       !                                   ! -------------------------------------------------- 
    358       zu_frc(:,:) = 0._wp 
    359       zv_frc(:,:) = 0._wp 
    360       ! 
    361       DO jk = 1, jpkm1 
    362          zu_frc(:,:) = zu_frc(:,:) + e3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 
    363          zv_frc(:,:) = zv_frc(:,:) + e3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)          
    364       END DO 
    365       ! 
    366       zu_frc(:,:) = zu_frc(:,:) * r1_hu_n(:,:) 
    367       zv_frc(:,:) = zv_frc(:,:) * r1_hv_n(:,:) 
    368       ! 
    369       ! 
    370       !                                   !* baroclinic momentum trend (remove the vertical mean trend) 
    371       DO jk = 1, jpkm1                    ! ----------------------------------------------------------- 
    372          DO jj = 2, jpjm1 
    373             DO ji = fs_2, fs_jpim1   ! vector opt. 
    374                ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk) 
    375                va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk) 
    376             END DO 
    377          END DO 
     227      !                                   !=  zu_frc =  1/H e3*d/dt(Ua)  =!  (Vertical mean of Ua, the 3D trends) 
     228      !                                   !  ---------------------------  ! 
     229      zu_frc(:,:) = SUM( e3u_n(:,:,:) * ua(:,:,:) * umask(:,:,:) , DIM=3 ) * r1_hu_n(:,:) 
     230      zv_frc(:,:) = SUM( e3v_n(:,:,:) * va(:,:,:) * vmask(:,:,:) , DIM=3 ) * r1_hv_n(:,:) 
     231      ! 
     232      ! 
     233      !                                   !=  Ua => baroclinic trend  =!   (remove its vertical mean) 
     234      DO jk = 1, jpkm1                    !  ------------------------  ! 
     235         ua(:,:,jk) = ( ua(:,:,jk) - zu_frc(:,:) ) * umask(:,:,jk) 
     236         va(:,:,jk) = ( va(:,:,jk) - zv_frc(:,:) ) * vmask(:,:,jk) 
    378237      END DO 
    379238       
     
    381240!!gm  Is it correct to do so ?   I think so... 
    382241       
    383        
    384       !                                   !* barotropic Coriolis trends (vorticity scheme dependent) 
    385       !                                   ! -------------------------------------------------------- 
    386       ! 
    387       zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:)        ! now fluxes  
    388       zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:) 
    389       ! 
    390       SELECT CASE( nvor_scheme ) 
    391       CASE( np_ENT )                ! enstrophy conserving scheme (f-point) 
    392          DO jj = 2, jpjm1 
    393             DO ji = 2, jpim1   ! vector opt. 
    394                zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * r1_hu_n(ji,jj)                    & 
    395                   &               * (  e1e2t(ji+1,jj)*ht_n(ji+1,jj)*ff_t(ji+1,jj) * ( vn_b(ji+1,jj) + vn_b(ji+1,jj-1) )   & 
    396                   &                  + e1e2t(ji  ,jj)*ht_n(ji  ,jj)*ff_t(ji  ,jj) * ( vn_b(ji  ,jj) + vn_b(ji  ,jj-1) )   ) 
    397                   ! 
    398                zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * r1_hv_n(ji,jj)                    & 
    399                   &               * (  e1e2t(ji,jj+1)*ht_n(ji,jj+1)*ff_t(ji,jj+1) * ( un_b(ji,jj+1) + un_b(ji-1,jj+1) )   &  
    400                   &                  + e1e2t(ji,jj  )*ht_n(ji,jj  )*ff_t(ji,jj  ) * ( un_b(ji,jj  ) + un_b(ji-1,jj  ) )   )  
    401             END DO   
    402          END DO   
    403          !          
    404       CASE( np_ENE , np_MIX )        ! energy conserving scheme (t-point) ENE or MIX 
    405          DO jj = 2, jpjm1 
    406             DO ji = fs_2, fs_jpim1   ! vector opt. 
    407                zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj) 
    408                zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj) 
    409                zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj) 
    410                zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj) 
    411                ! energy conserving formulation for planetary vorticity term 
    412                zu_trd(ji,jj) =   r1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
    413                zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
    414             END DO 
    415          END DO 
    416          ! 
    417       CASE( np_ENS )                ! enstrophy conserving scheme (f-point) 
    418          DO jj = 2, jpjm1 
    419             DO ji = fs_2, fs_jpim1   ! vector opt. 
    420                zy1 =   r1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) & 
    421                  &            + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj) 
    422                zx1 = - r1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) & 
    423                  &            + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj) 
    424                zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) ) 
    425                zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) ) 
    426             END DO 
    427          END DO 
    428          ! 
    429       CASE( np_EET , np_EEN )      ! energy & enstrophy scheme (using e3t or e3f)          
    430          DO jj = 2, jpjm1 
    431             DO ji = fs_2, fs_jpim1   ! vector opt. 
    432                zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) & 
    433                 &                                         + ftnw(ji+1,jj) * zwy(ji+1,jj  ) & 
    434                 &                                         + ftse(ji,jj  ) * zwy(ji  ,jj-1) & 
    435                 &                                         + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 
    436                zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 
    437                 &                                         + ftse(ji,jj+1) * zwx(ji  ,jj+1) & 
    438                 &                                         + ftnw(ji,jj  ) * zwx(ji-1,jj  ) & 
    439                 &                                         + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) 
    440             END DO 
    441          END DO 
    442          ! 
    443       END SELECT 
    444       ! 
    445       !                                   !* Right-Hand-Side of the barotropic momentum equation 
    446       !                                   ! ---------------------------------------------------- 
    447       IF( .NOT.ln_linssh ) THEN                 ! Variable volume : remove surface pressure gradient 
    448          IF( ln_wd_il ) THEN                        ! Calculating and applying W/D gravity filters 
     242      !                                   !=  remove 2D Coriolis and pressure gradient trends  =! 
     243      !                                   !  -------------------------------------------------  ! 
     244      ! 
     245      IF( kt == nit000 .OR. .NOT. ln_linssh )   CALL dyn_cor_2D_init   ! Set zwz, the barotropic Coriolis force coefficient 
     246      !       ! recompute zwz = f/depth  at every time step for (.NOT.ln_linssh) as the water colomn height changes 
     247      ! 
     248      !                                         !* 2D Coriolis trends 
     249      zhU(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:)        ! now fluxes  
     250      zhV(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:)        ! NB: FULL domain : put a value in last row and column 
     251      ! 
     252      CALL dyn_cor_2d( hu_n, hv_n, un_b, vn_b, zhU, zhV,  &   ! <<== in 
     253         &                               zu_trd, zv_trd   )   ! ==>> out 
     254      ! 
     255      IF( .NOT.ln_linssh ) THEN                 !* surface pressure gradient   (variable volume only) 
     256         ! 
     257         IF( ln_wd_il ) THEN                       ! W/D : limiter applied to spgspg 
     258            CALL wad_spg( sshn, zcpx, zcpy )          ! Calculating W/D gravity filters, zcpx and zcpy 
    449259            DO jj = 2, jpjm1 
    450                DO ji = 2, jpim1  
    451                   ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji+1,jj) ) >                & 
    452                      &      MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            & 
    453                      &      MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji+1,jj) + ht_0(ji+1,jj) )  & 
    454                      &                                                         > rn_wdmin1 + rn_wdmin2 
    455                   ll_tmp2 = ( ABS( sshn(ji+1,jj)            -  sshn(ji  ,jj))  > 1.E-12 ).AND.( & 
    456                      &      MAX(   sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                & 
    457                      &      MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
    458                   IF(ll_tmp1) THEN 
    459                      zcpx(ji,jj) = 1.0_wp 
    460                   ELSEIF(ll_tmp2) THEN 
    461                      ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
    462                      zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 
    463                                  &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) ) 
    464                      zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 
    465                   ELSE 
    466                      zcpx(ji,jj) = 0._wp 
    467                   ENDIF 
    468                   ! 
    469                   ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji,jj+1) ) >                & 
    470                      &      MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            & 
    471                      &      MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji,jj+1) + ht_0(ji,jj+1) )  & 
    472                      &                                                       > rn_wdmin1 + rn_wdmin2 
    473                   ll_tmp2 = ( ABS( sshn(ji,jj)              -  sshn(ji,jj+1))  > 1.E-12 ).AND.( & 
    474                      &      MAX(   sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                & 
    475                      &      MAX(  -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
    476    
    477                   IF(ll_tmp1) THEN 
    478                      zcpy(ji,jj) = 1.0_wp 
    479                   ELSE IF(ll_tmp2) THEN 
    480                      ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
    481                      zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 
    482                         &             / (sshn(ji,jj+1) - sshn(ji,jj  )) ) 
    483                      zcpy(ji,jj) = MAX(  0._wp , MIN( zcpy(ji,jj) , 1.0_wp )  ) 
    484                   ELSE 
    485                      zcpy(ji,jj) = 0._wp 
    486                   ENDIF 
    487                END DO 
    488             END DO 
    489             ! 
    490             DO jj = 2, jpjm1 
    491                DO ji = 2, jpim1 
     260               DO ji = 2, jpim1                ! SPG with the application of W/D gravity filters 
    492261                  zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   & 
    493262                     &                          * r1_e1u(ji,jj) * zcpx(ji,jj)  * wdrampu(ji,jj)  !jth 
     
    496265               END DO 
    497266            END DO 
    498             ! 
    499          ELSE 
    500             ! 
     267         ELSE                                      ! now suface pressure gradient 
    501268            DO jj = 2, jpjm1 
    502269               DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    516283      END DO  
    517284      ! 
    518       !                                         ! Add bottom stress contribution from baroclinic velocities:       
    519       IF (ln_bt_fw) THEN 
    520          DO jj = 2, jpjm1                           
    521             DO ji = fs_2, fs_jpim1   ! vector opt. 
    522                ikbu = mbku(ji,jj)        
    523                ikbv = mbkv(ji,jj)     
    524                zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities 
    525                zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj) 
    526             END DO 
    527          END DO 
    528       ELSE 
    529          DO jj = 2, jpjm1 
    530             DO ji = fs_2, fs_jpim1   ! vector opt. 
    531                ikbu = mbku(ji,jj)        
    532                ikbv = mbkv(ji,jj)     
    533                zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities 
    534                zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj) 
    535             END DO 
    536          END DO 
    537       ENDIF 
    538       ! 
    539       ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 
    540       IF( ln_wd_il ) THEN 
    541          zztmp = -1._wp / rdtbt 
    542          DO jj = 2, jpjm1 
    543             DO ji = fs_2, fs_jpim1   ! vector opt. 
    544                zu_frc(ji,jj) = zu_frc(ji,jj) + &  
    545                & MAX(r1_hu_n(ji,jj) * r1_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ), zztmp ) * zwx(ji,jj) *  wdrampu(ji,jj) 
    546                zv_frc(ji,jj) = zv_frc(ji,jj) + &  
    547                & MAX(r1_hv_n(ji,jj) * r1_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ), zztmp ) * zwy(ji,jj) *  wdrampv(ji,jj) 
    548             END DO 
    549          END DO 
    550       ELSE 
    551          DO jj = 2, jpjm1 
    552             DO ji = fs_2, fs_jpim1   ! vector opt. 
    553                zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * r1_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zwx(ji,jj) 
    554                zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * r1_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zwy(ji,jj) 
    555             END DO 
    556          END DO 
    557       END IF 
    558       ! 
    559       IF( ln_isfcav ) THEN       ! Add TOP stress contribution from baroclinic velocities:       
    560          IF( ln_bt_fw ) THEN 
    561             DO jj = 2, jpjm1 
     285      !                                   !=  Add bottom stress contribution from baroclinic velocities  =! 
     286      !                                   !  -----------------------------------------------------------  ! 
     287      CALL dyn_drg_init( zu_frc, zv_frc,  zCdU_u, zCdU_v )      ! also provide the barotropic drag coefficients 
     288      ! 
     289      !                                   !=  Add atmospheric pressure forcing  =! 
     290      !                                   !  ----------------------------------  ! 
     291      IF( ln_apr_dyn ) THEN 
     292         IF( ln_bt_fw ) THEN                          ! FORWARD integration: use kt+1/2 pressure (NOW+1/2) 
     293            DO jj = 2, jpjm1               
    562294               DO ji = fs_2, fs_jpim1   ! vector opt. 
    563                   iktu = miku(ji,jj) 
    564                   iktv = mikv(ji,jj) 
    565                   zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities 
    566                   zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj) 
    567                END DO 
    568             END DO 
    569          ELSE 
    570             DO jj = 2, jpjm1 
     295                  zu_frc(ji,jj) = zu_frc(ji,jj) + grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj) 
     296                  zv_frc(ji,jj) = zv_frc(ji,jj) + grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj) 
     297               END DO 
     298            END DO 
     299         ELSE                                         ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW) 
     300            zztmp = grav * r1_2 
     301            DO jj = 2, jpjm1               
    571302               DO ji = fs_2, fs_jpim1   ! vector opt. 
    572                   iktu = miku(ji,jj) 
    573                   iktv = mikv(ji,jj) 
    574                   zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities 
    575                   zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj) 
    576                END DO 
    577             END DO 
    578          ENDIF 
    579          ! 
    580          ! Note that the "unclipped" top friction parameter is used even with explicit drag 
    581          DO jj = 2, jpjm1               
    582             DO ji = fs_2, fs_jpim1   ! vector opt. 
    583                zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * r1_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zwx(ji,jj) 
    584                zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * r1_2 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zwy(ji,jj) 
    585             END DO 
    586          END DO 
    587       ENDIF 
    588       !        
     303                  zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)  & 
     304                       &                                   + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj) 
     305                  zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj)  & 
     306                       &                                   + ssh_ibb(ji  ,jj+1) - ssh_ibb(ji,jj)  ) * r1_e2v(ji,jj) 
     307               END DO 
     308            END DO 
     309         ENDIF  
     310      ENDIF 
     311      ! 
     312      !                                   !=  Add atmospheric pressure forcing  =! 
     313      !                                   !  ----------------------------------  ! 
    589314      IF( ln_bt_fw ) THEN                        ! Add wind forcing 
    590315         DO jj = 2, jpjm1 
     
    604329      ENDIF   
    605330      ! 
    606       IF( ln_apr_dyn ) THEN                     ! Add atm pressure forcing 
    607          IF( ln_bt_fw ) THEN 
    608             DO jj = 2, jpjm1               
    609                DO ji = fs_2, fs_jpim1   ! vector opt. 
    610                   zu_spg =  grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj) 
    611                   zv_spg =  grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj) 
    612                   zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 
    613                   zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg 
    614                END DO 
    615             END DO 
    616          ELSE 
    617             zztmp = grav * r1_2 
    618             DO jj = 2, jpjm1               
    619                DO ji = fs_2, fs_jpim1   ! vector opt. 
    620                   zu_spg = zztmp * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)    & 
    621                       &             + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj) 
    622                   zv_spg = zztmp * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj)    & 
    623                       &             + ssh_ibb(ji  ,jj+1) - ssh_ibb(ji,jj)  ) * r1_e2v(ji,jj) 
    624                   zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 
    625                   zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg 
    626                END DO 
    627             END DO 
    628          ENDIF  
    629       ENDIF 
    630       !                                   !* Right-Hand-Side of the barotropic ssh equation 
    631       !                                   ! ----------------------------------------------- 
    632       !                                         ! Surface net water flux and rivers 
    633       IF (ln_bt_fw) THEN 
    634          zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) 
    635       ELSE 
     331      !              !----------------! 
     332      !              !==  sssh_frc  ==!   Right-Hand-Side of the barotropic ssh equation   (over the FULL domain) 
     333      !              !----------------! 
     334      !                                   !=  Net water flux forcing applied to a water column  =! 
     335      !                                   ! ---------------------------------------------------  ! 
     336      IF (ln_bt_fw) THEN                          ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2) 
     337         zssh_frc(:,:) = r1_rau0 * ( emp(:,:)             - rnf(:,:)              + fwfisf(:,:)                  ) 
     338      ELSE                                        ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW) 
    636339         zztmp = r1_rau0 * r1_2 
    637          zssh_frc(:,:) = zztmp * (  emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)   & 
    638                 &                 + fwfisf(:,:) + fwfisf_b(:,:)                     ) 
    639       ENDIF 
    640       ! 
    641       IF( ln_sdw ) THEN                         ! Stokes drift divergence added if necessary 
     340         zssh_frc(:,:) = zztmp * (  emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) + fwfisf(:,:) + fwfisf_b(:,:)  ) 
     341      ENDIF 
     342      !                                   !=  Add Stokes drift divergence  =!   (if exist) 
     343      IF( ln_sdw ) THEN                   !  -----------------------------  ! 
    642344         zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:) 
    643345      ENDIF 
    644346      ! 
    645347#if defined key_asminc 
    646       !                                         ! Include the IAU weighted SSH increment 
     348      !                                   !=  Add the IAU weighted SSH increment  =! 
     349      !                                   !  ------------------------------------  ! 
    647350      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 
    648351         zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:) 
    649352      ENDIF 
    650353#endif 
    651       !                                   !* Fill boundary data arrays for AGRIF 
     354      !                                   != Fill boundary data arrays for AGRIF 
    652355      !                                   ! ------------------------------------ 
    653356#if defined key_agrif 
     
    671374         vb_e   (:,:) = 0._wp 
    672375      ENDIF 
    673  
     376      ! 
     377      IF( ln_linssh ) THEN    ! mid-step ocean depth is fixed (hup2_e=hu_n=hu_0) 
     378         zhup2_e(:,:) = hu_n(:,:) 
     379         zhvp2_e(:,:) = hv_n(:,:) 
     380         zhtp2_e(:,:) = ht_n(:,:) 
     381      ENDIF 
    674382      ! 
    675383      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                     
     
    693401      ENDIF 
    694402      ! 
    695       ! 
    696       ! 
    697403      ! Initialize sums: 
    698404      ua_b  (:,:) = 0._wp       ! After barotropic velocities (or transport if flux form)           
     
    714420         ! 
    715421         l_full_nf_update = jn == icycle   ! false: disable full North fold update (performances) for jn = 1 to icycle-1 
    716          !                                                !  ------------------ 
    717          !                                                !* Update the forcing (BDY and tides) 
    718          !                                                !  ------------------ 
    719          ! Update only tidal forcing at open boundaries 
    720          IF( ln_bdy      .AND. ln_tide )   CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 ) 
    721          IF( ln_tide_pot .AND. ln_tide )   CALL upd_tide     ( kt, kit=jn, time_offset= noffset   ) 
    722          ! 
    723          ! Set extrapolation coefficients for predictor step: 
     422         ! 
     423         !                    !==  Update the forcing ==! (BDY and tides) 
     424         ! 
     425         IF( ln_bdy      .AND. ln_tide )   CALL bdy_dta_tides( kt, kit=jn, kt_offset= noffset+1 ) 
     426         IF( ln_tide_pot .AND. ln_tide )   CALL upd_tide     ( kt, kit=jn, kt_offset= noffset   ) 
     427         ! 
     428         !                    !==  extrapolation at mid-step  ==!   (jn+1/2) 
     429         ! 
     430         !                       !* Set extrapolation coefficients for predictor step: 
    724431         IF ((jn<3).AND.ll_init) THEN      ! Forward            
    725432           za1 = 1._wp                                           
     
    731438           za3 =  0.281105_wp              ! za3 = bet 
    732439         ENDIF 
    733  
    734          ! Extrapolate barotropic velocities at step jit+0.5: 
     440         ! 
     441         !                       !* Extrapolate barotropic velocities at mid-step (jn+1/2) 
     442         !--        m+1/2               m                m-1           m-2       --! 
     443         !--       u      = (3/2+beta) u   -(1/2+2beta) u      + beta u          --! 
     444         !-------------------------------------------------------------------------! 
    735445         ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:) 
    736446         va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) 
     
    739449            !                                             !  ------------------ 
    740450            ! Extrapolate Sea Level at step jit+0.5: 
     451            !--         m+1/2                 m                  m-1             m-2       --! 
     452            !--      ssh      = (3/2+beta) ssh   -(1/2+2beta) ssh      + beta ssh          --! 
     453            !--------------------------------------------------------------------------------! 
    741454            zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 
    742455             
    743             ! set wetting & drying mask at tracer points for this barotropic sub-step  
    744             IF ( ln_wd_dl ) THEN  
    745                ! 
    746                IF ( ln_wd_dl_rmp ) THEN  
    747                   DO jj = 1, jpj                                  
    748                      DO ji = 1, jpi   ! vector opt.   
    749                         IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN  
    750 !                        IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  rn_wdmin2 ) THEN  
    751                            ztwdmask(ji,jj) = 1._wp 
    752                         ELSE IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN 
    753                            ztwdmask(ji,jj) = (tanh(50._wp*( ( zsshp2_e(ji,jj) + ht_0(ji,jj) -  rn_wdmin1 )*r_rn_wdmin1)) )  
    754                         ELSE  
    755                            ztwdmask(ji,jj) = 0._wp 
    756                         END IF 
    757                      END DO 
    758                   END DO 
    759                ELSE 
    760                   DO jj = 1, jpj                                  
    761                      DO ji = 1, jpi   ! vector opt.   
    762                         IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN  
    763                            ztwdmask(ji,jj) = 1._wp 
    764                         ELSE  
    765                            ztwdmask(ji,jj) = 0._wp 
    766                         ENDIF 
    767                      END DO 
    768                   END DO 
    769                ENDIF  
    770                ! 
    771             ENDIF  
     456            ! set wetting & drying mask at tracer points for this barotropic mid-step 
     457            IF( ln_wd_dl )   CALL wad_tmsk( zsshp2_e, ztwdmask ) 
    772458            ! 
    773             DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points 
    774                DO ji = 2, fs_jpim1   ! Vector opt. 
    775                   zwx(ji,jj) = r1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj)     & 
    776                      &              * ( e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  & 
    777                      &              +   e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) 
    778                   zwy(ji,jj) = r1_2 * ssvmask(ji,jj)  * r1_e1e2v(ji,jj)     & 
    779                      &              * ( e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )  & 
    780                      &              +   e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) 
    781                END DO 
    782             END DO 
    783             CALL lbc_lnk_multi( 'dynspg_ts', zwx, 'U', 1._wp, zwy, 'V', 1._wp ) 
     459            !                          ! ocean t-depth at mid-step 
     460            zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:) 
    784461            ! 
    785             zhup2_e(:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points 
    786             zhvp2_e(:,:) = hv_0(:,:) + zwy(:,:) 
    787             zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:) 
    788          ELSE 
    789             zhup2_e(:,:) = hu_n(:,:) 
    790             zhvp2_e(:,:) = hv_n(:,:) 
    791             zhtp2_e(:,:) = ht_n(:,:) 
    792          ENDIF 
    793          !                                                !* after ssh 
    794          !                                                !  ----------- 
    795          ! 
    796          ! Enforce volume conservation at open boundaries: 
     462            !                          ! ocean u- and v-depth at mid-step   (separate DO-loops remove the need of a lbc_lnk) 
     463            DO jj = 1, jpj 
     464               DO ji = 1, jpim1   ! not jpi-column 
     465                  zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj)                        & 
     466                       &                              * (  e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  & 
     467                       &                                 + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj)  ) * ssumask(ji,jj) 
     468               END DO 
     469            END DO 
     470            DO jj = 1, jpjm1        ! not jpj-row 
     471               DO ji = 1, jpi 
     472                  zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj)                        & 
     473                       &                              * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )  & 
     474                       &                                 + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1)  ) * ssvmask(ji,jj) 
     475               END DO 
     476            END DO 
     477            ! 
     478         ENDIF 
     479         ! 
     480         !                    !==  after SSH  ==!   (jn+1) 
     481         ! 
     482         !                             ! update (ua_e,va_e) to enforce volume conservation at open boundaries 
     483         !                             ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d 
    797484         IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e ) 
    798485         ! 
    799          zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)         ! fluxes at jn+0.5 
    800          zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 
     486         !                             ! resulting flux at mid-step (not over the full domain) 
     487         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 
     488         zhV(1:jpi  ,1:jpjm1) = e1v(1:jpi  ,1:jpjm1) * va_e(1:jpi  ,1:jpjm1) * zhvp2_e(1:jpi  ,1:jpjm1)   ! not jpj-row 
    801489         ! 
    802490#if defined key_agrif 
     
    805493            IF((nbondi == -1).OR.(nbondi == 2)) THEN 
    806494               DO jj = 1, jpj 
    807                   zwx(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * e2u(2:nbghostcells+1,jj) 
    808                   zwy(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * e1v(2:nbghostcells+1,jj) 
     495                  zhU(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * e2u(2:nbghostcells+1,jj) 
     496                  zhV(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * e1v(2:nbghostcells+1,jj) 
    809497               END DO 
    810498            ENDIF 
    811499            IF((nbondi ==  1).OR.(nbondi == 2)) THEN 
    812500               DO jj=1,jpj 
    813                   zwx(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj) 
    814                   zwy(nlci-nbghostcells  :nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * e1v(nlci-nbghostcells  :nlci-1,jj) 
     501                  zhU(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj) 
     502                  zhV(nlci-nbghostcells  :nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * e1v(nlci-nbghostcells  :nlci-1,jj) 
    815503               END DO 
    816504            ENDIF 
    817505            IF((nbondj == -1).OR.(nbondj == 2)) THEN 
    818506               DO ji=1,jpi 
    819                   zwy(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1) 
    820                   zwx(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * e2u(ji,2:nbghostcells+1) 
     507                  zhV(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1) 
     508                  zhU(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * e2u(ji,2:nbghostcells+1) 
    821509               END DO 
    822510            ENDIF 
    823511            IF((nbondj ==  1).OR.(nbondj == 2)) THEN 
    824512               DO ji=1,jpi 
    825                   zwy(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2) 
    826                   zwx(ji,nlcj-nbghostcells  :nlcj-1) = ubdy_n(ji,1:nbghostcells) * e2u(ji,nlcj-nbghostcells  :nlcj-1) 
     513                  zhV(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2) 
     514                  zhU(ji,nlcj-nbghostcells  :nlcj-1) = ubdy_n(ji,1:nbghostcells) * e2u(ji,nlcj-nbghostcells  :nlcj-1) 
    827515               END DO 
    828516            ENDIF 
    829517         ENDIF 
    830518#endif 
    831          IF( ln_wd_il )   CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 
    832  
    833          IF ( ln_wd_dl ) THEN  
    834             ! 
    835             ! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells  
    836             ! 
    837             DO jj = 1, jpjm1                                  
    838                DO ji = 1, jpim1    
    839                   IF ( zwx(ji,jj) > 0.0 ) THEN  
    840                      zuwdmask(ji, jj) = ztwdmask(ji  ,jj)  
    841                   ELSE  
    842                      zuwdmask(ji, jj) = ztwdmask(ji+1,jj)   
    843                   END IF  
    844                   zwx(ji, jj) = zuwdmask(ji,jj)*zwx(ji, jj) 
    845                   un_e(ji,jj) = zuwdmask(ji,jj)*un_e(ji,jj) 
    846  
    847                   IF ( zwy(ji,jj) > 0.0 ) THEN 
    848                      zvwdmask(ji, jj) = ztwdmask(ji, jj  ) 
    849                   ELSE  
    850                      zvwdmask(ji, jj) = ztwdmask(ji, jj+1)   
    851                   END IF  
    852                   zwy(ji, jj) = zvwdmask(ji,jj)*zwy(ji,jj)  
    853                   vn_e(ji,jj) = zvwdmask(ji,jj)*vn_e(ji,jj) 
    854                END DO 
    855             END DO 
     519         IF( ln_wd_il )   CALL wad_lmt_bt(zhU, zhV, sshn_e, zssh_frc, rdtbt)    !!gm wad_lmt_bt use of lbc_lnk on zhU, zhV 
     520 
     521         IF( ln_wd_dl ) THEN           ! un_e and vn_e are set to zero at faces where  
     522            !                          ! the direction of the flow is from dry cells 
     523            CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask )   ! not jpi colomn for U, not jpj row for V 
    856524            ! 
    857525         ENDIF     
    858           
    859          ! Sum over sub-time-steps to compute advective velocities 
    860          za2 = wgtbtp2(jn) 
    861          un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:) 
    862          vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:) 
    863           
    864          ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc = True)  
     526         ! 
     527         ! 
     528         !     Compute Sea Level at step jit+1 
     529         !--           m+1        m                               m+1/2          --! 
     530         !--        ssh    =  ssh   - delta_t' * [ frc + div( flux      ) ]      --! 
     531         !-------------------------------------------------------------------------! 
     532         DO jj = 2, jpjm1        ! INNER domain                              
     533            DO ji = 2, jpim1 
     534               zhdiv = (   zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1)   ) * r1_e1e2t(ji,jj) 
     535               ssha_e(ji,jj) = (  sshn_e(ji,jj) - rdtbt * ( zssh_frc(ji,jj) + zhdiv )  ) * ssmask(ji,jj) 
     536            END DO 
     537         END DO 
     538         ! 
     539         CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp,  zhU, 'U', -1._wp,  zhV, 'V', -1._wp ) 
     540         ! 
     541         !                             ! Sum over sub-time-steps to compute advective velocities 
     542         za2 = wgtbtp2(jn)             ! zhU, zhV hold fluxes extrapolated at jn+0.5 
     543         un_adv(:,:) = un_adv(:,:) + za2 * zhU(:,:) * r1_e2u(:,:) 
     544         vn_adv(:,:) = vn_adv(:,:) + za2 * zhV(:,:) * r1_e1v(:,:) 
     545         ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True)  
    865546         IF ( ln_wd_dl_bc ) THEN 
    866             zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:) 
    867             zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:) 
    868          END IF  
    869  
    870          ! Set next sea level: 
    871          DO jj = 2, jpjm1                                  
    872             DO ji = fs_2, fs_jpim1   ! vector opt. 
    873                zhdiv(ji,jj) = (   zwx(ji,jj) - zwx(ji-1,jj)   & 
    874                   &             + zwy(ji,jj) - zwy(ji,jj-1)   ) * r1_e1e2t(ji,jj) 
    875             END DO 
    876          END DO 
    877          ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * ssmask(:,:) 
    878           
    879          CALL lbc_lnk( 'dynspg_ts', ssha_e, 'T',  1._wp ) 
    880  
     547            zuwdav2(1:jpim1,1:jpj  ) = zuwdav2(1:jpim1,1:jpj  ) + za2 * zuwdmask(1:jpim1,1:jpj  )   ! not jpi-column 
     548            zvwdav2(1:jpi  ,1:jpjm1) = zvwdav2(1:jpi  ,1:jpjm1) + za2 * zvwdmask(1:jpi  ,1:jpjm1)   ! not jpj-row 
     549         END IF 
     550         ! 
    881551         ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T) 
    882552         IF( ln_bdy )   CALL bdy_ssh( ssha_e ) 
     
    887557         ! Sea Surface Height at u-,v-points (vvl case only) 
    888558         IF( .NOT.ln_linssh ) THEN                                 
    889             DO jj = 2, jpjm1 
     559            DO jj = 2, jpjm1   ! INNER domain, will be extended to whole domain later 
    890560               DO ji = 2, jpim1      ! NO Vector Opt. 
    891561                  zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj)    & 
     
    897567               END DO 
    898568            END DO 
    899             CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) 
    900569         ENDIF    
    901          !                                  
    902          ! Half-step back interpolation of SSH for surface pressure computation: 
    903          !---------------------------------------------------------------------- 
    904          IF ((jn==1).AND.ll_init) THEN 
    905            za0=1._wp                        ! Forward-backward 
    906            za1=0._wp                            
    907            za2=0._wp 
    908            za3=0._wp 
    909          ELSEIF ((jn==2).AND.ll_init) THEN  ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12 
    910            za0= 1.0833333333333_wp          ! za0 = 1-gam-eps 
    911            za1=-0.1666666666666_wp          ! za1 = gam 
    912            za2= 0.0833333333333_wp          ! za2 = eps 
    913            za3= 0._wp               
    914          ELSE                               ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880  
    915             IF (rn_bt_alpha==0._wp) THEN 
    916                za0=0.614_wp                 ! za0 = 1/2 +   gam + 2*eps 
    917                za1=0.285_wp                 ! za1 = 1/2 - 2*gam - 3*eps 
    918                za2=0.088_wp                 ! za2 = gam 
    919                za3=0.013_wp                 ! za3 = eps 
    920             ELSE 
    921                zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha 
    922                zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha 
    923                za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon 
    924                za1 = 1._wp - za0 - zgamma - zepsilon 
    925                za2 = zgamma 
    926                za3 = zepsilon 
    927             ENDIF  
    928          ENDIF 
    929          ! 
     570         !          
     571         ! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2 
     572         !--            m+1/2           m+1              m               m-1              m-2     --! 
     573         !--        ssh'    =  za0 * ssh     +  za1 * ssh   +  za2 * ssh      +  za3 * ssh        --! 
     574         !------------------------------------------------------------------------------------------! 
     575         CALL ts_bck_interp( jn, ll_init, za0, za1, za2, za3 )   ! coeficients of the interpolation 
    930576         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:)   & 
    931577            &          + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:) 
    932           
    933          IF( ln_wd_il ) THEN                   ! Calculating and applying W/D gravity filters 
    934            DO jj = 2, jpjm1 
    935               DO ji = 2, jpim1  
    936                 ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji+1,jj) ) >                & 
    937                      &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji+1,jj) ) .AND.            & 
    938                      &    MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) & 
    939                      &                                                             > rn_wdmin1 + rn_wdmin2 
    940                 ll_tmp2 = (ABS(zsshp2_e(ji,jj)               - zsshp2_e(ji+1,jj))  > 1.E-12 ).AND.( & 
    941                      &    MAX( zsshp2_e(ji,jj)               , zsshp2_e(ji+1,jj) ) >                & 
    942                      &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
    943     
    944                 IF(ll_tmp1) THEN 
    945                   zcpx(ji,jj) = 1.0_wp 
    946                 ELSE IF(ll_tmp2) THEN 
    947                   ! no worries about  zsshp2_e(ji+1,jj) - zsshp2_e(ji  ,jj) = 0, it won't happen ! here 
    948                   zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) +     ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 
    949                               &    / (zsshp2_e(ji+1,jj) - zsshp2_e(ji  ,jj)) ) 
    950                 ELSE 
    951                   zcpx(ji,jj) = 0._wp 
    952                 ENDIF 
    953                 ! 
    954                 ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                & 
    955                      &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji,jj+1) ) .AND.            & 
    956                      &    MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) & 
    957                      &                                                             > rn_wdmin1 + rn_wdmin2 
    958                 ll_tmp2 = (ABS(zsshp2_e(ji,jj)               - zsshp2_e(ji,jj+1))  > 1.E-12 ).AND.( & 
    959                      &    MAX( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                & 
    960                      &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
    961     
    962                 IF(ll_tmp1) THEN 
    963                   zcpy(ji,jj) = 1.0_wp 
    964                 ELSEIF(ll_tmp2) THEN 
    965                   ! no worries about  zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  ) = 0, it won't happen ! here 
    966                   zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) +     ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 
    967                               &    / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  )) ) 
    968                 ELSE 
    969                   zcpy(ji,jj) = 0._wp 
    970                 ENDIF 
    971               END DO 
    972            END DO 
    973          ENDIF 
    974          ! 
    975          ! Compute associated depths at U and V points: 
    976          IF( .NOT.ln_linssh  .AND. .NOT.ln_dynadv_vec ) THEN     !* Vector form 
    977             !                                         
    978             DO jj = 2, jpjm1                             
    979                DO ji = 2, jpim1 
    980                   zx1 = r1_2 * ssumask(ji  ,jj) *  r1_e1e2u(ji  ,jj)    & 
    981                      &      * ( e1e2t(ji  ,jj  ) * zsshp2_e(ji  ,jj)    & 
    982                      &      +   e1e2t(ji+1,jj  ) * zsshp2_e(ji+1,jj  ) ) 
    983                   zy1 = r1_2 * ssvmask(ji  ,jj) *  r1_e1e2v(ji  ,jj  )  & 
    984                      &       * ( e1e2t(ji ,jj  ) * zsshp2_e(ji  ,jj  )  & 
    985                      &       +   e1e2t(ji ,jj+1) * zsshp2_e(ji  ,jj+1) ) 
    986                   zhust_e(ji,jj) = hu_0(ji,jj) + zx1  
    987                   zhvst_e(ji,jj) = hv_0(ji,jj) + zy1 
    988                END DO 
    989             END DO 
    990             ! 
     578         ! 
     579         !                             ! Surface pressure gradient 
     580         zldg = ( 1._wp - rn_scal_load ) * grav    ! local factor 
     581         DO jj = 2, jpjm1                             
     582            DO ji = 2, jpim1 
     583               zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 
     584               zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 
     585            END DO 
     586         END DO 
     587         IF( ln_wd_il ) THEN        ! W/D : gravity filters applied on pressure gradient 
     588            CALL wad_spg( zsshp2_e, zcpx, zcpy )   ! Calculating W/D gravity filters 
     589            zu_spg(2:jpim1,2:jpjm1) = zu_spg(2:jpim1,2:jpjm1) * zcpx(2:jpim1,2:jpjm1) 
     590            zv_spg(2:jpim1,2:jpjm1) = zv_spg(2:jpim1,2:jpjm1) * zcpy(2:jpim1,2:jpjm1) 
    991591         ENDIF 
    992592         ! 
     
    994594         ! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated 
    995595         ! at each time step. We however keep them constant here for optimization. 
    996          ! Recall that zwx and zwy arrays hold fluxes at this stage: 
    997          ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)   ! fluxes at jn+0.5 
    998          ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 
    999          ! 
    1000          SELECT CASE( nvor_scheme ) 
    1001          CASE( np_ENT )             ! energy conserving scheme (t-point) 
    1002          DO jj = 2, jpjm1 
    1003             DO ji = 2, jpim1   ! vector opt. 
    1004  
    1005                z1_hu = ssumask(ji,jj) / ( zhup2_e(ji,jj) + 1._wp - ssumask(ji,jj) ) 
    1006                z1_hv = ssvmask(ji,jj) / ( zhvp2_e(ji,jj) + 1._wp - ssvmask(ji,jj) ) 
    1007              
    1008                zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu                   & 
    1009                   &               * (  e1e2t(ji+1,jj)*zhtp2_e(ji+1,jj)*ff_t(ji+1,jj) * ( va_e(ji+1,jj) + va_e(ji+1,jj-1) )   & 
    1010                   &                  + e1e2t(ji  ,jj)*zhtp2_e(ji  ,jj)*ff_t(ji  ,jj) * ( va_e(ji  ,jj) + va_e(ji  ,jj-1) )   ) 
    1011                   ! 
    1012                zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv                    & 
    1013                   &               * (  e1e2t(ji,jj+1)*zhtp2_e(ji,jj+1)*ff_t(ji,jj+1) * ( ua_e(ji,jj+1) + ua_e(ji-1,jj+1) )   &  
    1014                   &                  + e1e2t(ji,jj  )*zhtp2_e(ji,jj  )*ff_t(ji,jj  ) * ( ua_e(ji,jj  ) + ua_e(ji-1,jj  ) )   )  
    1015             END DO   
    1016          END DO   
    1017          !          
    1018          CASE( np_ENE, np_MIX )     ! energy conserving scheme (f-point) 
    1019             DO jj = 2, jpjm1 
    1020                DO ji = fs_2, fs_jpim1   ! vector opt. 
    1021                   zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj) 
    1022                   zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj) 
    1023                   zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj) 
    1024                   zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj) 
    1025                   zu_trd(ji,jj) = r1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
    1026                   zv_trd(ji,jj) =-r1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
    1027                END DO 
    1028             END DO 
    1029             ! 
    1030          CASE( np_ENS )             ! enstrophy conserving scheme (f-point) 
    1031             DO jj = 2, jpjm1 
    1032                DO ji = fs_2, fs_jpim1   ! vector opt. 
    1033                   zy1 =   r1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) & 
    1034                    &             + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj) 
    1035                   zx1 = - r1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) & 
    1036                    &             + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj) 
    1037                   zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) ) 
    1038                   zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) ) 
    1039                END DO 
    1040             END DO 
    1041             ! 
    1042          CASE( np_EET , np_EEN )   ! energy & enstrophy scheme (using e3t or e3f) 
    1043             DO jj = 2, jpjm1 
    1044                DO ji = fs_2, fs_jpim1   ! vector opt. 
    1045                   zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  )  & 
    1046                      &                                       + ftnw(ji+1,jj) * zwy(ji+1,jj  )  & 
    1047                      &                                       + ftse(ji,jj  ) * zwy(ji  ,jj-1)  &  
    1048                      &                                       + ftsw(ji+1,jj) * zwy(ji+1,jj-1)  ) 
    1049                   zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1)  &  
    1050                      &                                       + ftse(ji,jj+1) * zwx(ji  ,jj+1)  & 
    1051                      &                                       + ftnw(ji,jj  ) * zwx(ji-1,jj  )  &  
    1052                      &                                       + ftne(ji,jj  ) * zwx(ji  ,jj  )  ) 
    1053                END DO 
    1054             END DO 
    1055             !  
    1056          END SELECT 
     596         ! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated) 
     597         CALL dyn_cor_2d( zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV,    zu_trd, zv_trd   ) 
    1057598         ! 
    1058599         ! Add tidal astronomical forcing if defined 
     
    1060601            DO jj = 2, jpjm1 
    1061602               DO ji = fs_2, fs_jpim1   ! vector opt. 
    1062                   zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 
    1063                   zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 
    1064                   zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg 
    1065                   zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg 
     603                  zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 
     604                  zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 
    1066605               END DO 
    1067606            END DO 
     
    1077616               END DO 
    1078617            END DO 
    1079          ENDIF  
    1080          ! 
    1081          ! Surface pressure trend: 
    1082          IF( ln_wd_il ) THEN 
    1083            DO jj = 2, jpjm1 
    1084               DO ji = 2, jpim1  
    1085                  ! Add surface pressure gradient 
    1086                  zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 
    1087                  zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 
    1088                  zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg * zcpx(ji,jj)  
    1089                  zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg * zcpy(ji,jj) 
    1090               END DO 
    1091            END DO 
    1092          ELSE 
    1093            DO jj = 2, jpjm1 
    1094               DO ji = fs_2, fs_jpim1   ! vector opt. 
    1095                  ! Add surface pressure gradient 
    1096                  zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 
    1097                  zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 
    1098                  zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg 
    1099                  zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg 
    1100               END DO 
    1101            END DO 
    1102          END IF 
    1103  
     618         ENDIF 
    1104619         ! 
    1105620         ! Set next velocities: 
     621         !     Compute barotropic speeds at step jit+1    (h : total height of the water colomn) 
     622         !--                              VECTOR FORM 
     623         !--   m+1                 m               /                                                       m+1/2           \    --! 
     624         !--  u     =             u   + delta_t' * \         (1-r)*g * grad_x( ssh') -         f * k vect u      +     frc /    --! 
     625         !--                                                                                                                    --! 
     626         !--                             FLUX FORM                                                                              --! 
     627         !--  m+1   __1__  /  m    m               /  m+1/2                             m+1/2              m+1/2    n      \ \  --! 
     628         !-- u    =   m+1 |  h  * u   + delta_t' * \ h     * (1-r)*g * grad_x( ssh') - h     * f * k vect u      + h * frc /  | --! 
     629         !--         h     \                                                                                                 /  --! 
     630         !------------------------------------------------------------------------------------------------------------------------! 
    1106631         IF( ln_dynadv_vec .OR. ln_linssh ) THEN      !* Vector form 
    1107632            DO jj = 2, jpjm1 
    1108633               DO ji = fs_2, fs_jpim1   ! vector opt. 
    1109634                  ua_e(ji,jj) = (                                 un_e(ji,jj)   &  
    1110                             &     + rdtbt * (                      zwx(ji,jj)   & 
     635                            &     + rdtbt * (                   zu_spg(ji,jj)   & 
    1111636                            &                                 + zu_trd(ji,jj)   & 
    1112637                            &                                 + zu_frc(ji,jj) ) &  
     
    1114639 
    1115640                  va_e(ji,jj) = (                                 vn_e(ji,jj)   & 
    1116                             &     + rdtbt * (                      zwy(ji,jj)   & 
     641                            &     + rdtbt * (                   zv_spg(ji,jj)   & 
    1117642                            &                                 + zv_trd(ji,jj)   & 
    1118643                            &                                 + zv_frc(ji,jj) ) & 
    1119644                            &   ) * ssvmask(ji,jj) 
    1120   
    1121645               END DO 
    1122646            END DO 
     
    1124648         ELSE                           !* Flux form 
    1125649            DO jj = 2, jpjm1 
    1126                DO ji = fs_2, fs_jpim1   ! vector opt. 
    1127  
    1128                   zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 
    1129                   zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 
    1130  
    1131                   zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj)) 
    1132                   zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj)) 
    1133  
    1134                   ua_e(ji,jj) = (                hu_e(ji,jj)  *   un_e(ji,jj)   &  
    1135                             &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   &  
    1136                             &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   & 
    1137                             &               +    hu_n(ji,jj)  * zu_frc(ji,jj) ) & 
    1138                             &   ) * zhura 
    1139  
    1140                   va_e(ji,jj) = (                hv_e(ji,jj)  *   vn_e(ji,jj)   & 
    1141                             &     + rdtbt * ( zhvst_e(ji,jj)  *    zwy(ji,jj)   & 
    1142                             &               + zhvp2_e(ji,jj)  * zv_trd(ji,jj)   & 
    1143                             &               +    hv_n(ji,jj)  * zv_frc(ji,jj) ) & 
    1144                             &   ) * zhvra 
     650               DO ji = 2, jpim1 
     651                  !                    ! hu_e, hv_e hold depth at jn,  zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2 
     652                  !                    ! backward interpolated depth used in spg terms at jn+1/2 
     653                  zhu_bck = hu_0(ji,jj) + r1_2*r1_e1e2u(ji,jj) * (  e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)    & 
     654                       &                                          + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj)  ) * ssumask(ji,jj) 
     655                  zhv_bck = hv_0(ji,jj) + r1_2*r1_e1e2v(ji,jj) * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )    & 
     656                       &                                          + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1)  ) * ssvmask(ji,jj) 
     657                  !                    ! inverse depth at jn+1 
     658                  z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) ) 
     659                  z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) ) 
     660                  ! 
     661                  ua_e(ji,jj) = (               hu_e  (ji,jj) *   un_e (ji,jj)      &  
     662                       &            + rdtbt * (  zhu_bck        * zu_spg (ji,jj)  &   ! 
     663                       &                       + zhup2_e(ji,jj) * zu_trd (ji,jj)  &   ! 
     664                       &                       +  hu_n  (ji,jj) * zu_frc (ji,jj)  )   ) * z1_hu 
     665                  ! 
     666                  va_e(ji,jj) = (               hv_e  (ji,jj) *   vn_e (ji,jj)      & 
     667                       &            + rdtbt * (  zhv_bck        * zv_spg (ji,jj)  &   ! 
     668                       &                       + zhvp2_e(ji,jj) * zv_trd (ji,jj)  &   ! 
     669                       &                       +  hv_n  (ji,jj) * zv_frc (ji,jj)  )   ) * z1_hv 
    1145670               END DO 
    1146671            END DO 
     
    1155680            END DO 
    1156681         ENDIF 
    1157  
    1158           
    1159          IF( .NOT.ln_linssh ) THEN                     !* Update ocean depth (variable volume case only) 
    1160             hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 
    1161             hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 
    1162             hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) ) 
    1163             hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) ) 
    1164             ! 
    1165          ENDIF 
    1166          !                                             !* domain lateral boundary 
    1167          CALL lbc_lnk_multi( 'dynspg_ts', ua_e, 'U', -1._wp, va_e , 'V', -1._wp ) 
     682        
     683         IF( .NOT.ln_linssh ) THEN   !* Update ocean depth (variable volume case only) 
     684            hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1) 
     685            hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) ) 
     686            hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1) 
     687            hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) ) 
     688            CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp  & 
     689                 &                         , hu_e , 'U', -1._wp, hv_e , 'V', -1._wp  & 
     690                 &                         , hur_e, 'U', -1._wp, hvr_e, 'V', -1._wp  ) 
     691         ELSE 
     692            CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp  ) 
     693         ENDIF 
     694         ! 
    1168695         ! 
    1169696         !                                                 ! open boundaries 
     
    1213740      ! Set advection velocity correction: 
    1214741      IF (ln_bt_fw) THEN 
    1215          zwx(:,:) = un_adv(:,:) 
    1216          zwy(:,:) = vn_adv(:,:) 
    1217742         IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN 
    1218             un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) ) 
    1219             vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) ) 
    1220             ! 
    1221             ! Update corrective fluxes for next time step: 
    1222             un_bf(:,:)  = atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:)) 
    1223             vn_bf(:,:)  = atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:)) 
     743            DO jj = 1, jpj 
     744               DO ji = 1, jpi 
     745                  zun_save = un_adv(ji,jj) 
     746                  zvn_save = vn_adv(ji,jj) 
     747                  !                          ! apply the previously computed correction  
     748                  un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - atfp * un_bf(ji,jj) ) 
     749                  vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - atfp * vn_bf(ji,jj) ) 
     750                  !                          ! Update corrective fluxes for next time step 
     751                  un_bf(ji,jj)  = atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) ) 
     752                  vn_bf(ji,jj)  = atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) ) 
     753                  !                          ! Save integrated transport for next computation 
     754                  ub2_b(ji,jj) = zun_save 
     755                  vb2_b(ji,jj) = zvn_save 
     756               END DO 
     757            END DO 
    1224758         ELSE 
    1225             un_bf(:,:) = 0._wp 
    1226             vn_bf(:,:) = 0._wp  
    1227          END IF          
    1228          ! Save integrated transport for next computation 
    1229          ub2_b(:,:) = zwx(:,:) 
    1230          vb2_b(:,:) = zwy(:,:) 
     759            un_bf(:,:) = 0._wp            ! corrective fluxes for next time step set to zero 
     760            vn_bf(:,:) = 0._wp 
     761            ub2_b(:,:) = un_adv(:,:)      ! Save integrated transport for next computation 
     762            vb2_b(:,:) = vn_adv(:,:) 
     763         END IF 
    1231764      ENDIF 
    1232765 
     
    14731006      REAL(wp) ::   zxr2, zyr2, zcmax   ! local scalar 
    14741007      REAL(wp), DIMENSION(jpi,jpj) ::   zcu 
     1008      INTEGER  :: inum 
    14751009      !!---------------------------------------------------------------------- 
    14761010      ! 
     
    15791113   END SUBROUTINE dyn_spg_ts_init 
    15801114 
     1115    
     1116   SUBROUTINE dyn_cor_2d_init 
     1117      !!--------------------------------------------------------------------- 
     1118      !!                   ***  ROUTINE dyn_cor_2d_init  *** 
     1119      !! 
     1120      !! ** Purpose : Set time splitting options 
     1121      !! Set arrays to remove/compute coriolis trend. 
     1122      !! Do it once during initialization if volume is fixed, else at each long time step. 
     1123      !! Note that these arrays are also used during barotropic loop. These are however frozen 
     1124      !! although they should be updated in the variable volume case. Not a big approximation. 
     1125      !! To remove this approximation, copy lines below inside barotropic loop 
     1126      !! and update depths at T-F points (ht and zhf resp.) at each barotropic time step 
     1127      !! 
     1128      !! Compute zwz = f / ( height of the water colomn ) 
     1129      !!---------------------------------------------------------------------- 
     1130      INTEGER  ::   ji ,jj, jk              ! dummy loop indices 
     1131      REAL(wp) ::   z1_ht 
     1132      REAL(wp), DIMENSION(jpi,jpj) :: zhf 
     1133      !!---------------------------------------------------------------------- 
     1134      ! 
     1135      SELECT CASE( nvor_scheme ) 
     1136      CASE( np_EEN )                != EEN scheme using e3f (energy & enstrophy scheme) 
     1137         SELECT CASE( nn_een_e3f )              !* ff_f/e3 at F-point 
     1138         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4) 
     1139            DO jj = 1, jpjm1 
     1140               DO ji = 1, jpim1 
     1141                  zwz(ji,jj) =   ( ht_n(ji  ,jj+1) + ht_n(ji+1,jj+1) +                    & 
     1142                       &             ht_n(ji  ,jj  ) + ht_n(ji+1,jj  )   ) * 0.25_wp   
     1143                  IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 
     1144               END DO 
     1145            END DO 
     1146         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask) 
     1147            DO jj = 1, jpjm1 
     1148               DO ji = 1, jpim1 
     1149                  zwz(ji,jj) =             (  ht_n  (ji  ,jj+1) + ht_n  (ji+1,jj+1)      & 
     1150                       &                      + ht_n  (ji  ,jj  ) + ht_n  (ji+1,jj  )  )   & 
     1151                       &       / ( MAX( 1._wp,  ssmask(ji  ,jj+1) + ssmask(ji+1,jj+1)      & 
     1152                       &                      + ssmask(ji  ,jj  ) + ssmask(ji+1,jj  )  )   ) 
     1153                  IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 
     1154               END DO 
     1155            END DO 
     1156         END SELECT 
     1157         CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp ) 
     1158         ! 
     1159         ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 
     1160         DO jj = 2, jpj 
     1161            DO ji = 2, jpi 
     1162               ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1) 
     1163               ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  ) 
     1164               ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1) 
     1165               ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  ) 
     1166            END DO 
     1167         END DO 
     1168         ! 
     1169      CASE( np_EET )                  != EEN scheme using e3t (energy conserving scheme) 
     1170         ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 
     1171         DO jj = 2, jpj 
     1172            DO ji = 2, jpi 
     1173               z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) ) 
     1174               ftne(ji,jj) = ( ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) ) * z1_ht 
     1175               ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) ) * z1_ht 
     1176               ftse(ji,jj) = ( ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht 
     1177               ftsw(ji,jj) = ( ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) ) * z1_ht 
     1178            END DO 
     1179         END DO 
     1180         ! 
     1181      CASE( np_ENE, np_ENS , np_MIX )  != all other schemes (ENE, ENS, MIX) except ENT ! 
     1182         ! 
     1183         zwz(:,:) = 0._wp 
     1184         zhf(:,:) = 0._wp 
     1185          
     1186         !!gm  assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed  
     1187!!gm    A priori a better value should be something like : 
     1188!!gm          zhf(i,j) = masked sum of  ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)  
     1189!!gm                     divided by the sum of the corresponding mask  
     1190!!gm  
     1191!!             
     1192         IF( .NOT.ln_sco ) THEN 
     1193   
     1194   !!gm  agree the JC comment  : this should be done in a much clear way 
     1195   
     1196   ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 
     1197   !     Set it to zero for the time being  
     1198   !              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level 
     1199   !              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth 
     1200   !              ENDIF 
     1201   !              zhf(:,:) = gdepw_0(:,:,jk+1) 
     1202            ! 
     1203         ELSE 
     1204            ! 
     1205            !zhf(:,:) = hbatf(:,:) 
     1206            DO jj = 1, jpjm1 
     1207               DO ji = 1, jpim1 
     1208                  zhf(ji,jj) =    (   ht_0  (ji,jj  ) + ht_0  (ji+1,jj  )          & 
     1209                       &            + ht_0  (ji,jj+1) + ht_0  (ji+1,jj+1)   )      & 
     1210                       &     / MAX(   ssmask(ji,jj  ) + ssmask(ji+1,jj  )          & 
     1211                       &            + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp  ) 
     1212               END DO 
     1213            END DO 
     1214         ENDIF 
     1215         ! 
     1216         DO jj = 1, jpjm1 
     1217            zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 
     1218         END DO 
     1219         ! 
     1220         DO jk = 1, jpkm1 
     1221            DO jj = 1, jpjm1 
     1222               zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 
     1223            END DO 
     1224         END DO 
     1225         CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp ) 
     1226         ! JC: TBC. hf should be greater than 0  
     1227         DO jj = 1, jpj 
     1228            DO ji = 1, jpi 
     1229               IF( zhf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zhf(ji,jj) 
     1230            END DO 
     1231         END DO 
     1232         zwz(:,:) = ff_f(:,:) * zwz(:,:) 
     1233      END SELECT 
     1234       
     1235   END SUBROUTINE dyn_cor_2d_init 
     1236 
     1237 
     1238 
     1239   SUBROUTINE dyn_cor_2d( hu_n, hv_n, un_b, vn_b, zhU, zhV,    zu_trd, zv_trd   ) 
     1240      !!--------------------------------------------------------------------- 
     1241      !!                   ***  ROUTINE dyn_cor_2d  *** 
     1242      !! 
     1243      !! ** Purpose : Compute u and v coriolis trends 
     1244      !!---------------------------------------------------------------------- 
     1245      INTEGER  ::   ji ,jj                             ! dummy loop indices 
     1246      REAL(wp) ::   zx1, zx2, zy1, zy2, z1_hu, z1_hv   !   -      - 
     1247      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: hu_n, hv_n, un_b, vn_b, zhU, zhV 
     1248      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: zu_trd, zv_trd 
     1249      !!---------------------------------------------------------------------- 
     1250      SELECT CASE( nvor_scheme ) 
     1251      CASE( np_ENT )                ! enstrophy conserving scheme (f-point) 
     1252         DO jj = 2, jpjm1 
     1253            DO ji = 2, jpim1 
     1254               z1_hu = ssumask(ji,jj) / ( hu_n(ji,jj) + 1._wp - ssumask(ji,jj) ) 
     1255               z1_hv = ssvmask(ji,jj) / ( hv_n(ji,jj) + 1._wp - ssvmask(ji,jj) ) 
     1256               zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu                    & 
     1257                  &               * (  e1e2t(ji+1,jj)*ht_n(ji+1,jj)*ff_t(ji+1,jj) * ( vn_b(ji+1,jj) + vn_b(ji+1,jj-1) )   & 
     1258                  &                  + e1e2t(ji  ,jj)*ht_n(ji  ,jj)*ff_t(ji  ,jj) * ( vn_b(ji  ,jj) + vn_b(ji  ,jj-1) )   ) 
     1259                  ! 
     1260               zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv                    & 
     1261                  &               * (  e1e2t(ji,jj+1)*ht_n(ji,jj+1)*ff_t(ji,jj+1) * ( un_b(ji,jj+1) + un_b(ji-1,jj+1) )   &  
     1262                  &                  + e1e2t(ji,jj  )*ht_n(ji,jj  )*ff_t(ji,jj  ) * ( un_b(ji,jj  ) + un_b(ji-1,jj  ) )   )  
     1263            END DO   
     1264         END DO   
     1265         !          
     1266      CASE( np_ENE , np_MIX )        ! energy conserving scheme (t-point) ENE or MIX 
     1267         DO jj = 2, jpjm1 
     1268            DO ji = fs_2, fs_jpim1   ! vector opt. 
     1269               zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj) 
     1270               zy2 = ( zhV(ji,jj  ) + zhV(ji+1,jj  ) ) * r1_e1u(ji,jj) 
     1271               zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj) 
     1272               zx2 = ( zhU(ji  ,jj) + zhU(ji  ,jj+1) ) * r1_e2v(ji,jj) 
     1273               ! energy conserving formulation for planetary vorticity term 
     1274               zu_trd(ji,jj) =   r1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
     1275               zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
     1276            END DO 
     1277         END DO 
     1278         ! 
     1279      CASE( np_ENS )                ! enstrophy conserving scheme (f-point) 
     1280         DO jj = 2, jpjm1 
     1281            DO ji = fs_2, fs_jpim1   ! vector opt. 
     1282               zy1 =   r1_8 * ( zhV(ji  ,jj-1) + zhV(ji+1,jj-1) & 
     1283                 &            + zhV(ji  ,jj  ) + zhV(ji+1,jj  ) ) * r1_e1u(ji,jj) 
     1284               zx1 = - r1_8 * ( zhU(ji-1,jj  ) + zhU(ji-1,jj+1) & 
     1285                 &            + zhU(ji  ,jj  ) + zhU(ji  ,jj+1) ) * r1_e2v(ji,jj) 
     1286               zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) ) 
     1287               zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) ) 
     1288            END DO 
     1289         END DO 
     1290         ! 
     1291      CASE( np_EET , np_EEN )      ! energy & enstrophy scheme (using e3t or e3f)          
     1292         DO jj = 2, jpjm1 
     1293            DO ji = fs_2, fs_jpim1   ! vector opt. 
     1294               zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zhV(ji  ,jj  ) & 
     1295                &                                         + ftnw(ji+1,jj) * zhV(ji+1,jj  ) & 
     1296                &                                         + ftse(ji,jj  ) * zhV(ji  ,jj-1) & 
     1297                &                                         + ftsw(ji+1,jj) * zhV(ji+1,jj-1) ) 
     1298               zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zhU(ji-1,jj+1) & 
     1299                &                                         + ftse(ji,jj+1) * zhU(ji  ,jj+1) & 
     1300                &                                         + ftnw(ji,jj  ) * zhU(ji-1,jj  ) & 
     1301                &                                         + ftne(ji,jj  ) * zhU(ji  ,jj  ) ) 
     1302            END DO 
     1303         END DO 
     1304         ! 
     1305      END SELECT 
     1306      ! 
     1307   END SUBROUTINE dyn_cor_2D 
     1308 
     1309 
     1310   SUBROUTINE wad_tmsk( pssh, ptmsk ) 
     1311      !!---------------------------------------------------------------------- 
     1312      !!                  ***  ROUTINE wad_lmt  *** 
     1313      !!                     
     1314      !! ** Purpose :   set wetting & drying mask at tracer points  
     1315      !!              for the current barotropic sub-step  
     1316      !! 
     1317      !! ** Method  :   ???  
     1318      !! 
     1319      !! ** Action  :  ptmsk : wetting & drying t-mask 
     1320      !!---------------------------------------------------------------------- 
     1321      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pssh    ! 
     1322      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   ptmsk   ! 
     1323      ! 
     1324      INTEGER  ::   ji, jj   ! dummy loop indices 
     1325      !!---------------------------------------------------------------------- 
     1326      ! 
     1327      IF( ln_wd_dl_rmp ) THEN      
     1328         DO jj = 1, jpj 
     1329            DO ji = 1, jpi                     
     1330               IF    ( pssh(ji,jj) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN  
     1331                  !           IF    ( pssh(ji,jj) + ht_0(ji,jj) >          rn_wdmin2 ) THEN  
     1332                  ptmsk(ji,jj) = 1._wp 
     1333               ELSEIF( pssh(ji,jj) + ht_0(ji,jj) >          rn_wdmin1 ) THEN 
     1334                  ptmsk(ji,jj) = TANH( 50._wp*( ( pssh(ji,jj) + ht_0(ji,jj) -  rn_wdmin1 )*r_rn_wdmin1) ) 
     1335               ELSE  
     1336                  ptmsk(ji,jj) = 0._wp 
     1337               ENDIF 
     1338            END DO 
     1339         END DO 
     1340      ELSE   
     1341         DO jj = 1, jpj 
     1342            DO ji = 1, jpi                               
     1343               IF ( pssh(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN   ;   ptmsk(ji,jj) = 1._wp 
     1344               ELSE                                                 ;   ptmsk(ji,jj) = 0._wp 
     1345               ENDIF 
     1346            END DO 
     1347         END DO 
     1348      ENDIF 
     1349      ! 
     1350   END SUBROUTINE wad_tmsk 
     1351 
     1352 
     1353   SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk ) 
     1354      !!---------------------------------------------------------------------- 
     1355      !!                  ***  ROUTINE wad_lmt  *** 
     1356      !!                     
     1357      !! ** Purpose :   set wetting & drying mask at tracer points  
     1358      !!              for the current barotropic sub-step  
     1359      !! 
     1360      !! ** Method  :   ???  
     1361      !! 
     1362      !! ** Action  :  ptmsk : wetting & drying t-mask 
     1363      !!---------------------------------------------------------------------- 
     1364      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pTmsk              ! W & D t-mask 
     1365      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   phU, phV, pu, pv   ! ocean velocities and transports 
     1366      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pUmsk, pVmsk       ! W & D u- and v-mask 
     1367      ! 
     1368      INTEGER  ::   ji, jj   ! dummy loop indices 
     1369      !!---------------------------------------------------------------------- 
     1370      ! 
     1371      DO jj = 1, jpj 
     1372         DO ji = 1, jpim1   ! not jpi-column 
     1373            IF ( phU(ji,jj) > 0._wp ) THEN   ;   pUmsk(ji,jj) = pTmsk(ji  ,jj)  
     1374            ELSE                             ;   pUmsk(ji,jj) = pTmsk(ji+1,jj)   
     1375            ENDIF 
     1376            phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj) 
     1377            pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj) 
     1378         END DO 
     1379      END DO 
     1380      ! 
     1381      DO jj = 1, jpjm1   ! not jpj-row 
     1382         DO ji = 1, jpi 
     1383            IF ( phV(ji,jj) > 0._wp ) THEN   ;   pVmsk(ji,jj) = pTmsk(ji,jj  ) 
     1384            ELSE                             ;   pVmsk(ji,jj) = pTmsk(ji,jj+1)   
     1385            ENDIF 
     1386            phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj)  
     1387            pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj) 
     1388         END DO 
     1389      END DO 
     1390      ! 
     1391   END SUBROUTINE wad_Umsk 
     1392 
     1393 
     1394   SUBROUTINE wad_spg( sshn, zcpx, zcpy ) 
     1395      !!--------------------------------------------------------------------- 
     1396      !!                   ***  ROUTINE  wad_sp  *** 
     1397      !! 
     1398      !! ** Purpose :  
     1399      !!---------------------------------------------------------------------- 
     1400      INTEGER  ::   ji ,jj               ! dummy loop indices 
     1401      LOGICAL  ::   ll_tmp1, ll_tmp2 
     1402      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: sshn 
     1403      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy 
     1404      !!---------------------------------------------------------------------- 
     1405      DO jj = 2, jpjm1 
     1406         DO ji = 2, jpim1  
     1407            ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji+1,jj) ) >                & 
     1408                 &      MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            & 
     1409                 &      MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji+1,jj) + ht_0(ji+1,jj) )  & 
     1410                 &                                                         > rn_wdmin1 + rn_wdmin2 
     1411            ll_tmp2 = ( ABS( sshn(ji+1,jj)            -  sshn(ji  ,jj))  > 1.E-12 ).AND.( & 
     1412                 &      MAX(   sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                & 
     1413                 &      MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
     1414            IF(ll_tmp1) THEN 
     1415               zcpx(ji,jj) = 1.0_wp 
     1416            ELSEIF(ll_tmp2) THEN 
     1417               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
     1418               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 
     1419                    &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) ) 
     1420               zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 
     1421            ELSE 
     1422               zcpx(ji,jj) = 0._wp 
     1423            ENDIF 
     1424            ! 
     1425            ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji,jj+1) ) >                & 
     1426                 &      MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            & 
     1427                 &      MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji,jj+1) + ht_0(ji,jj+1) )  & 
     1428                 &                                                       > rn_wdmin1 + rn_wdmin2 
     1429            ll_tmp2 = ( ABS( sshn(ji,jj)              -  sshn(ji,jj+1))  > 1.E-12 ).AND.( & 
     1430                 &      MAX(   sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                & 
     1431                 &      MAX(  -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
     1432             
     1433            IF(ll_tmp1) THEN 
     1434               zcpy(ji,jj) = 1.0_wp 
     1435            ELSE IF(ll_tmp2) THEN 
     1436               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
     1437               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 
     1438                    &             / (sshn(ji,jj+1) - sshn(ji,jj  )) ) 
     1439               zcpy(ji,jj) = MAX(  0._wp , MIN( zcpy(ji,jj) , 1.0_wp )  ) 
     1440            ELSE 
     1441               zcpy(ji,jj) = 0._wp 
     1442            ENDIF 
     1443         END DO 
     1444      END DO 
     1445             
     1446   END SUBROUTINE wad_spg 
     1447      
     1448 
     1449 
     1450   SUBROUTINE dyn_drg_init( pu_RHSi, pv_RHSi, pCdU_u, pCdU_v ) 
     1451      !!---------------------------------------------------------------------- 
     1452      !!                  ***  ROUTINE dyn_drg_init  *** 
     1453      !!                     
     1454      !! ** Purpose : - add the baroclinic top/bottom drag contribution to  
     1455      !!              the baroclinic part of the barotropic RHS 
     1456      !!              - compute the barotropic drag coefficients 
     1457      !! 
     1458      !! ** Method  :   computation done over the INNER domain only  
     1459      !!---------------------------------------------------------------------- 
     1460      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pu_RHSi, pv_RHSi   ! baroclinic part of the barotropic RHS 
     1461      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pCdU_u , pCdU_v    ! barotropic drag coefficients 
     1462      ! 
     1463      INTEGER  ::   ji, jj   ! dummy loop indices 
     1464      INTEGER  ::   ikbu, ikbv, iktu, iktv 
     1465      REAL(wp) ::   zztmp 
     1466      REAL(wp), DIMENSION(jpi,jpj) ::   zu_i, zv_i 
     1467      !!---------------------------------------------------------------------- 
     1468      ! 
     1469      !                    !==  Set the barotropic drag coef.  ==! 
     1470      ! 
     1471      IF( ln_isfcav ) THEN          ! top+bottom friction (ocean cavities) 
     1472          
     1473         DO jj = 2, jpjm1 
     1474            DO ji = 2, jpim1     ! INNER domain 
     1475               pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) 
     1476               pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) 
     1477            END DO 
     1478         END DO 
     1479      ELSE                          ! bottom friction only 
     1480         DO jj = 2, jpjm1 
     1481            DO ji = 2, jpim1  ! INNER domain 
     1482               pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) 
     1483               pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) 
     1484            END DO 
     1485         END DO 
     1486      ENDIF 
     1487      ! 
     1488      !                    !==  BOTTOM stress contribution from baroclinic velocities  ==! 
     1489      ! 
     1490      IF( ln_bt_fw ) THEN                 ! FORWARD integration: use NOW bottom baroclinic velocities 
     1491          
     1492         DO jj = 2, jpjm1 
     1493            DO ji = 2, jpim1  ! INNER domain 
     1494               ikbu = mbku(ji,jj)        
     1495               ikbv = mbkv(ji,jj)     
     1496               zu_i(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) 
     1497               zv_i(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj) 
     1498            END DO 
     1499         END DO 
     1500      ELSE                                ! CENTRED integration: use BEFORE bottom baroclinic velocities 
     1501          
     1502         DO jj = 2, jpjm1 
     1503            DO ji = 2, jpim1   ! INNER domain 
     1504               ikbu = mbku(ji,jj)        
     1505               ikbv = mbkv(ji,jj)     
     1506               zu_i(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) 
     1507               zv_i(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj) 
     1508            END DO 
     1509         END DO 
     1510      ENDIF 
     1511      ! 
     1512      IF( ln_wd_il ) THEN      ! W/D : use the "clipped" bottom friction   !!gm   explain WHY, please ! 
     1513         zztmp = -1._wp / rdtbt 
     1514         DO jj = 2, jpjm1 
     1515            DO ji = 2, jpim1    ! INNER domain 
     1516               pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) *  wdrampu(ji,jj) * MAX(                                 &  
     1517                    &                              r1_hu_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp  ) 
     1518               pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) *  wdrampv(ji,jj) * MAX(                                 &  
     1519                    &                              r1_hv_n(ji,jj) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp  ) 
     1520            END DO 
     1521         END DO 
     1522      ELSE                    ! use "unclipped" drag (even if explicit friction is used in 3D calculation) 
     1523          
     1524         DO jj = 2, jpjm1 
     1525            DO ji = 2, jpim1    ! INNER domain 
     1526               pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj) 
     1527               pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv_n(ji,jj) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zv_i(ji,jj) 
     1528            END DO 
     1529         END DO 
     1530      END IF 
     1531      ! 
     1532      !                    !==  TOP stress contribution from baroclinic velocities  ==!   (no W/D case) 
     1533      ! 
     1534      IF( ln_isfcav ) THEN 
     1535         ! 
     1536         IF( ln_bt_fw ) THEN                ! FORWARD integration: use NOW top baroclinic velocity 
     1537             
     1538            DO jj = 2, jpjm1 
     1539               DO ji = 2, jpim1   ! INNER domain 
     1540                  iktu = miku(ji,jj) 
     1541                  iktv = mikv(ji,jj) 
     1542                  zu_i(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) 
     1543                  zv_i(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj) 
     1544               END DO 
     1545            END DO 
     1546         ELSE                                ! CENTRED integration: use BEFORE top baroclinic velocity 
     1547             
     1548            DO jj = 2, jpjm1 
     1549               DO ji = 2, jpim1      ! INNER domain 
     1550                  iktu = miku(ji,jj) 
     1551                  iktv = mikv(ji,jj) 
     1552                  zu_i(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) 
     1553                  zv_i(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj) 
     1554               END DO 
     1555            END DO 
     1556         ENDIF 
     1557         ! 
     1558         !                    ! use "unclipped" top drag (even if explicit friction is used in 3D calculation) 
     1559          
     1560         DO jj = 2, jpjm1 
     1561            DO ji = 2, jpim1    ! INNER domain 
     1562               pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu_n(ji,jj) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj) 
     1563               pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv_n(ji,jj) * r1_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zv_i(ji,jj) 
     1564            END DO 
     1565         END DO 
     1566         ! 
     1567      ENDIF 
     1568      ! 
     1569   END SUBROUTINE dyn_drg_init 
     1570 
     1571   SUBROUTINE ts_bck_interp( jn, ll_init,       &   ! <== in 
     1572      &                      za0, za1, za2, za3 )   ! ==> out 
     1573      !!---------------------------------------------------------------------- 
     1574      INTEGER ,INTENT(in   ) ::   jn                   ! index of sub time step 
     1575      LOGICAL ,INTENT(in   ) ::   ll_init              ! 
     1576      REAL(wp),INTENT(  out) ::   za0, za1, za2, za3   ! Half-step back interpolation coefficient 
     1577      ! 
     1578      REAL(wp) ::   zepsilon, zgamma                   !   -      - 
     1579      !!---------------------------------------------------------------------- 
     1580      !                             ! set Half-step back interpolation coefficient 
     1581      IF    ( jn==1 .AND. ll_init ) THEN   !* Forward-backward 
     1582         za0 = 1._wp                         
     1583         za1 = 0._wp                            
     1584         za2 = 0._wp 
     1585         za3 = 0._wp 
     1586      ELSEIF( jn==2 .AND. ll_init ) THEN   !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12 
     1587         za0 = 1.0833333333333_wp                 ! za0 = 1-gam-eps 
     1588         za1 =-0.1666666666666_wp                 ! za1 = gam 
     1589         za2 = 0.0833333333333_wp                 ! za2 = eps 
     1590         za3 = 0._wp               
     1591      ELSE                                 !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880  
     1592         IF( rn_bt_alpha == 0._wp ) THEN      ! Time diffusion   
     1593            za0 = 0.614_wp                        ! za0 = 1/2 +   gam + 2*eps 
     1594            za1 = 0.285_wp                        ! za1 = 1/2 - 2*gam - 3*eps 
     1595            za2 = 0.088_wp                        ! za2 = gam 
     1596            za3 = 0.013_wp                        ! za3 = eps 
     1597         ELSE                                 ! no time diffusion 
     1598            zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha 
     1599            zgamma   = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha 
     1600            za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon 
     1601            za1 = 1._wp - za0 - zgamma - zepsilon 
     1602            za2 = zgamma 
     1603            za3 = zepsilon 
     1604         ENDIF  
     1605      ENDIF 
     1606   END SUBROUTINE ts_bck_interp 
     1607 
     1608 
    15811609   !!====================================================================== 
    15821610END MODULE dynspg_ts 
Note: See TracChangeset for help on using the changeset viewer.