Changeset 11234


Ignore:
Timestamp:
2019-07-09T18:09:16+02:00 (13 months ago)
Author:
girrmann
Message:

dev_r10984_HPC-13 : cleaning of dyn_spg_ts routine, does not change the results, see #2285

Location:
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdydyn3d.F90

    r11210 r11234  
    186186      igrd = 3                      ! Copying tangential velocity into bdy points 
    187187      IF( llrim0 ) THEN   ;   ibeg = 1                       ;   iend = idx%nblenrim0(igrd) 
    188       ELSE               ;   ibeg = idx%nblenrim0(igrd)+1   ;   iend = idx%nblenrim(igrd) 
     188      ELSE                ;   ibeg = idx%nblenrim0(igrd)+1   ;   iend = idx%nblenrim(igrd) 
    189189      ENDIF 
    190190      DO jb = ibeg, iend 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyvol.F90

    r11049 r11234  
    9999            ii = idx%nbi(jb,jgrd) 
    100100            ij = idx%nbj(jb,jgrd) 
    101             IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove 
     101            IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! sum : else halo couted twice 
    102102            zubtpecor = zubtpecor + idx%flagu(jb,jgrd) * pua2d(ii,ij) * e2u(ii,ij) * phu(ii,ij) * tmask_i(ii,ij) * tmask_i(ii+1,ij) 
    103103         END DO 
     
    106106            ii = idx%nbi(jb,jgrd) 
    107107            ij = idx%nbj(jb,jgrd) 
    108             IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove 
     108            IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! sum : else halo couted twice 
    109109            zubtpecor = zubtpecor + idx%flagv(jb,jgrd) * pva2d(ii,ij) * e1v(ii,ij) * phv(ii,ij) * tmask_i(ii,ij) * tmask_i(ii,ij+1) 
    110110         END DO 
     
    128128               ii = idx%nbi(jb,jgrd) 
    129129               ij = idx%nbj(jb,jgrd) 
    130                IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! sum : else halo couted twice 
     130               !IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove ? 
    131131               pua2d(ii,ij) = pua2d(ii,ij) - idx%flagu(jb,jgrd) * zubtpecor * tmask_i(ii,ij) * tmask_i(ii+1,ij) 
    132132         END DO 
     
    135135               ii = idx%nbi(jb,jgrd) 
    136136               ij = idx%nbj(jb,jgrd) 
    137                IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! sum : else halo couted twice 
     137               !IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove ? 
    138138               pva2d(ii,ij) = pva2d(ii,ij) - idx%flagv(jb,jgrd) * zubtpecor * tmask_i(ii,ij) * tmask_i(ii,ij+1) 
    139139         END DO 
     
    154154                  ii = idx%nbi(jb,jgrd) 
    155155                  ij = idx%nbj(jb,jgrd) 
    156                   IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove? 
     156                  IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE 
    157157                  ztranst = ztranst + idx%flagu(jb,jgrd) * pua2d(ii,ij) * e2u(ii,ij) * phu(ii,ij) * tmask_i(ii,ij) * tmask_i(ii+1,ij) 
    158158            END DO 
     
    161161                  ii = idx%nbi(jb,jgrd) 
    162162                  ij = idx%nbj(jb,jgrd) 
    163                   IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove? 
     163                  IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE 
    164164                  ztranst = ztranst + idx%flagv(jb,jgrd) * pva2d(ii,ij) * e1v(ii,ij) * phv(ii,ij) * tmask_i(ii,ij) * tmask_i(ii,ij+1) 
    165165            END DO 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DYN/dynspg_ts.F90

    r11223 r11234  
    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) ) 
     
    153153      INTEGER  ::   ikbv, iktv            !   -      - 
    154154      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  !   -      - 
     155      REAL(wp) ::   zx1, zx2, zhura        , z1_hu  !   -      - 
     156      REAL(wp) ::   zy1, zy2, zhvra        , z1_hv  !   -      - 
    157157      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 
     158      REAL(wp) ::   zmdi, zztmp, zldg      , z1_ht  !   -      - 
     159      REAL(wp) ::   zhu_bck, zhv_bck                !   -      - 
     160      REAL(wp) ::   zun_save, zvn_save              !   -      - 
     161      REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e 
     162      REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zu_spg, zssh_frc 
     163      REAL(wp), DIMENSION(jpi,jpj) :: zwy, zv_trd, zv_frc, zv_spg, zhdiv 
    162164      REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhust_e, zhtp2_e 
    163165      REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zhvst_e 
    164166      REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v   ! top/bottom stress at u- & v-points 
     167      REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV         ! fluxes 
    165168      ! 
    166169      REAL(wp) ::   zwdramp                     ! local scalar - only used if ln_wd_dl = .True.  
     
    182185      zwdramp = r_rn_wdmin1               ! simplest ramp  
    183186!     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  
     187      !                                         ! inverse of baroclinic time step  
     188      IF( kt == nit000 .AND. neuler == 0 ) THEN   ;   r1_2dt_b = 1._wp / (         rdt ) 
     189      ELSE                                        ;   r1_2dt_b = 1._wp / ( 2._wp * rdt ) 
     190      ENDIF 
    189191      ! 
    190192      ll_init     = ln_bt_av                    ! if no time averaging, then no specific restart  
     
    210212            ll_fw_start =.FALSE. 
    211213         ENDIF 
    212          ! 
    213          ! Set averaging weights and cycle length: 
     214         !                    ! Set averaging weights and cycle length: 
    214215         CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 
    215216         ! 
    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 
    342217      ENDIF 
    343218      ! 
     
    348223         CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 
    349224      ENDIF 
     225      ! 
    350226                           
    351227      ! ----------------------------------------------------------------------------- 
     
    354230      !       
    355231      ! 
    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 
     232      !                                   !=  zu_frc =  1/H e3*d/dt(Ua)  =!  (Vertical mean of Ua, the 3D trends) 
     233      !                                   !  ---------------------------  ! 
     234      zu_frc(:,:) = SUM( e3u_n(:,:,:) * ua(:,:,:) * umask(:,:,:) , DIM=3 ) * r1_hu_n(:,:) 
     235      zv_frc(:,:) = SUM( e3v_n(:,:,:) * va(:,:,:) * vmask(:,:,:) , DIM=3 ) * r1_hv_n(:,:) 
     236      ! 
     237      ! 
     238      !                                   !=  Ua => baroclinic trend  =!   (remove its vertical mean) 
     239      DO jk = 1, jpkm1                    !  ------------------------  ! 
     240         ua(:,:,jk) = ( ua(:,:,jk) - zu_frc(:,:) ) * umask(:,:,jk) 
     241         va(:,:,jk) = ( va(:,:,jk) - zv_frc(:,:) ) * vmask(:,:,jk) 
    378242      END DO 
    379243       
     
    381245!!gm  Is it correct to do so ?   I think so... 
    382246       
    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 
     247      !                                   !=  remove 2D Coriolis and pressure gradient trends  =! 
     248      !                                   !  -------------------------------------------------  ! 
     249      ! 
     250      IF( kt == nit000 .OR. .NOT. ln_linssh )   CALL dyn_cor_2D_init   ! Set zwz, the barotropic Coriolis force coefficient 
     251      !       ! recompute zwz = f/depth  at every time step for (.NOT.ln_linssh) as the water colomn height changes 
     252      ! 
     253      !                                         !* 2D Coriolis trends 
     254      zhU(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:)        ! now fluxes  
     255      zhV(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:)        ! NB: FULL domain : put a value in last row and column 
     256      ! 
     257      CALL dyn_cor_2d( ht_n, hu_n, hv_n, un_b, vn_b, zhU, zhV,  &   ! <<== in 
     258         &                                     zu_trd, zv_trd   )   ! ==>> out 
     259      ! 
     260      IF( .NOT.ln_linssh ) THEN                 !* surface pressure gradient   (variable volume only) 
     261         ! 
     262         IF( ln_wd_il ) THEN                       ! W/D : limiter applied to spgspg 
     263            CALL wad_spg( sshn, zcpx, zcpy )          ! Calculating W/D gravity filters, zcpx and zcpy 
    449264            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 
     265               DO ji = 2, jpim1                ! SPG with the application of W/D gravity filters 
    492266                  zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   & 
    493267                     &                          * r1_e1u(ji,jj) * zcpx(ji,jj)  * wdrampu(ji,jj)  !jth 
     
    496270               END DO 
    497271            END DO 
    498             ! 
    499          ELSE 
    500             ! 
     272         ELSE                                      ! now suface pressure gradient 
    501273            DO jj = 2, jpjm1 
    502274               DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    516288      END DO  
    517289      ! 
    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 
     290      !                                   !=  Add bottom stress contribution from baroclinic velocities  =! 
     291      !                                   !  -----------------------------------------------------------  ! 
     292      CALL drg_init( zu_frc, zv_frc,  zCdU_u, zCdU_v )      ! also provide the barotropic drag coefficients 
     293      ! 
     294      !                                   !=  Add atmospheric pressure forcing  =! 
     295      !                                   !  ----------------------------------  ! 
     296      IF( ln_apr_dyn ) THEN 
     297         IF( ln_bt_fw ) THEN                          ! FORWARD integration: use kt+1/2 pressure (NOW+1/2) 
     298            DO jj = 2, jpjm1               
    562299               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 
     300                  zu_frc(ji,jj) = zu_frc(ji,jj) + grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj) 
     301                  zv_frc(ji,jj) = zv_frc(ji,jj) + grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj) 
     302               END DO 
     303            END DO 
     304         ELSE                                         ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW) 
     305            zztmp = grav * r1_2 
     306            DO jj = 2, jpjm1               
    571307               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       !        
     308                  zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)  & 
     309                       &                                   + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj) 
     310                  zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj)  & 
     311                       &                                   + ssh_ibb(ji  ,jj+1) - ssh_ibb(ji,jj)  ) * r1_e2v(ji,jj) 
     312               END DO 
     313            END DO 
     314         ENDIF  
     315      ENDIF 
     316      ! 
     317      !                                   !=  Add atmospheric pressure forcing  =! 
     318      !                                   !  ----------------------------------  ! 
    589319      IF( ln_bt_fw ) THEN                        ! Add wind forcing 
    590320         DO jj = 2, jpjm1 
     
    604334      ENDIF   
    605335      ! 
    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 
     336      !              !----------------! 
     337      !              !==  sssh_frc  ==!   Right-Hand-Side of the barotropic ssh equation   (over the FULL domain) 
     338      !              !----------------! 
     339      !                                   !=  Net water flux forcing applied to a water column  =! 
     340      !                                   ! ---------------------------------------------------  ! 
     341      IF (ln_bt_fw) THEN                          ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2) 
     342         zssh_frc(:,:) = r1_rau0 * ( emp(:,:)             - rnf(:,:)              + fwfisf(:,:)                  ) 
     343      ELSE                                        ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW) 
    636344         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 
     345         zssh_frc(:,:) = zztmp * (  emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) + fwfisf(:,:) + fwfisf_b(:,:)  ) 
     346      ENDIF 
     347      !                                   !=  Add Stokes drift divergence  =!   (if exist) 
     348      IF( ln_sdw ) THEN                   !  -----------------------------  ! 
    642349         zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:) 
    643350      ENDIF 
    644351      ! 
    645352#if defined key_asminc 
    646       !                                         ! Include the IAU weighted SSH increment 
     353      !                                   !=  Add the IAU weighted SSH increment  =! 
     354      !                                   !  ------------------------------------  ! 
    647355      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 
    648356         zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:) 
    649357      ENDIF 
    650358#endif 
    651       !                                   !* Fill boundary data arrays for AGRIF 
     359      !                                   != Fill boundary data arrays for AGRIF 
    652360      !                                   ! ------------------------------------ 
    653361#if defined key_agrif 
     
    671379         vb_e   (:,:) = 0._wp 
    672380      ENDIF 
    673  
     381      ! 
     382      IF( ln_linssh ) THEN    ! mid-step ocean depth is fixed (hup2_e=hu_n=hu_0) 
     383         zhup2_e(:,:) = hu_n(:,:) 
     384         zhvp2_e(:,:) = hv_n(:,:) 
     385         zhtp2_e(:,:) = ht_n(:,:) 
     386      ENDIF 
    674387      ! 
    675388      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                     
     
    693406      ENDIF 
    694407      ! 
    695       ! 
    696       ! 
    697408      ! Initialize sums: 
    698409      ua_b  (:,:) = 0._wp       ! After barotropic velocities (or transport if flux form)           
     
    714425         ! 
    715426         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 
     427         ! 
     428         !                    !==  Update the forcing ==! (BDY and tides) 
     429         ! 
    720430         IF( ln_bdy      .AND. ln_tide )   CALL bdy_dta_tides( kt, kit=jn, kt_offset= noffset+1 ) 
    721431         IF( ln_tide_pot .AND. ln_tide )   CALL upd_tide     ( kt, kit=jn, kt_offset= noffset   ) 
    722432         ! 
    723          ! Set extrapolation coefficients for predictor step: 
     433         !                    !==  extrapolation at mid-step  ==!   (jn+1/2) 
     434         ! 
     435         !                       !* Set extrapolation coefficients for predictor step: 
    724436         IF ((jn<3).AND.ll_init) THEN      ! Forward            
    725437           za1 = 1._wp                                           
     
    731443           za3 =  0.281105_wp              ! za3 = bet 
    732444         ENDIF 
    733  
    734          ! Extrapolate barotropic velocities at step jit+0.5: 
     445         ! 
     446         !                       !* Extrapolate barotropic velocities at mid-step (jn+1/2) 
     447         !--        m+1/2               m                m-1           m-2       --! 
     448         !--       u      = (3/2+beta) u   -(1/2+2beta) u      + beta u          --! 
     449         !-------------------------------------------------------------------------! 
    735450         ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:) 
    736451         va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) 
     
    739454            !                                             !  ------------------ 
    740455            ! Extrapolate Sea Level at step jit+0.5: 
     456            !--         m+1/2                 m                  m-1             m-2       --! 
     457            !--      ssh      = (3/2+beta) ssh   -(1/2+2beta) ssh      + beta ssh          --! 
     458            !--------------------------------------------------------------------------------! 
    741459            zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 
    742460             
    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  
     461            ! set wetting & drying mask at tracer points for this barotropic mid-step 
     462            IF( ln_wd_dl )   CALL wad_tmsk( zsshp2_e, ztwdmask ) 
    772463            ! 
    773464            DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points 
     
    786477            zhvp2_e(:,:) = hv_0(:,:) + zwy(:,:) 
    787478            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: 
     479            ! 
     480         ENDIF 
     481         ! 
     482         !                    !==  after SSH  ==!   (jn+1) 
     483         ! 
     484         !                             ! update (ua_e,va_e) to enforce volume conservation at open boundaries 
     485         !                             ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d 
    797486         IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e ) 
    798487         ! 
    799          zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)         ! fluxes at jn+0.5 
    800          zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 
     488         zhU(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)         ! fluxes at jn+0.5 
     489         zhV(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 
    801490         ! 
    802491#if defined key_agrif 
     
    805494            IF((nbondi == -1).OR.(nbondi == 2)) THEN 
    806495               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) 
     496                  zhU(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * e2u(2:nbghostcells+1,jj) 
     497                  zhV(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * e1v(2:nbghostcells+1,jj) 
    809498               END DO 
    810499            ENDIF 
    811500            IF((nbondi ==  1).OR.(nbondi == 2)) THEN 
    812501               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) 
     502                  zhU(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj) 
     503                  zhV(nlci-nbghostcells  :nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * e1v(nlci-nbghostcells  :nlci-1,jj) 
    815504               END DO 
    816505            ENDIF 
    817506            IF((nbondj == -1).OR.(nbondj == 2)) THEN 
    818507               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) 
     508                  zhV(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1) 
     509                  zhU(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * e2u(ji,2:nbghostcells+1) 
    821510               END DO 
    822511            ENDIF 
    823512            IF((nbondj ==  1).OR.(nbondj == 2)) THEN 
    824513               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) 
     514                  zhV(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2) 
     515                  zhU(ji,nlcj-nbghostcells  :nlcj-1) = ubdy_n(ji,1:nbghostcells) * e2u(ji,nlcj-nbghostcells  :nlcj-1) 
    827516               END DO 
    828517            ENDIF 
    829518         ENDIF 
    830519#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 
     520         IF( ln_wd_il )   CALL wad_lmt_bt(zhU, zhV, sshn_e, zssh_frc, rdtbt) 
     521 
     522         IF( ln_wd_dl ) THEN           ! un_e and vn_e are set to zero at faces where  
     523            !                          ! the direction of the flow is from dry cells 
     524            CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask ) 
    856525            ! 
    857526         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)  
     527         ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True)  
    865528         IF ( ln_wd_dl_bc ) THEN 
    866529            zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:) 
    867530            zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:) 
    868531         END IF  
     532          
     533         ! Sum over sub-time-steps to compute advective velocities 
     534         !  
     535         za2 = wgtbtp2(jn) 
     536         un_adv(:,:) = un_adv(:,:) + za2 * zhU(:,:) * r1_e2u(:,:) 
     537         vn_adv(:,:) = vn_adv(:,:) + za2 * zhV(:,:) * r1_e1v(:,:) 
    869538 
    870539         ! Set next sea level: 
    871540         DO jj = 2, jpjm1                                  
    872541            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 
     542               zhdiv(ji,jj) = (   zhU(ji,jj) - zhU(ji-1,jj)   & 
     543                  &             + zhV(ji,jj) - zhV(ji,jj-1)   ) * r1_e1e2t(ji,jj) 
     544            END DO 
     545         END DO 
     546         !     Compute Sea Level at step jit+1 
     547         !--           m+1        m                               m+1/2          --! 
     548         !--        ssh    =  ssh   - delta_t' * [ frc + div( flux      ) ]      --! 
     549         !-------------------------------------------------------------------------! 
    877550         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * ssmask(:,:) 
    878551          
     
    899572            CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) 
    900573         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          ! 
     574         !          
     575         ! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2 
     576         !--            m+1/2           m+1              m               m-1              m-2     --! 
     577         !--        ssh'    =  za0 * ssh     +  za1 * ssh   +  za2 * ssh      +  za3 * ssh        --! 
     578         !------------------------------------------------------------------------------------------! 
     579         CALL ts_bck_interp( jn, ll_init, za0, za1, za2, za3 )   ! coeficients of the interpolation 
    930580         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:)   & 
    931581            &          + 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             ! 
     582         ! 
     583         !                             ! Surface pressure gradient 
     584         zldg = ( 1._wp - rn_scal_load ) * grav    ! local factor 
     585         DO jj = 2, jpjm1                             
     586            DO ji = 2, jpim1 
     587               zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 
     588               zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 
     589            END DO 
     590         END DO 
     591         IF( ln_wd_il ) THEN        ! W/D : gravity filters applied on pressure gradient 
     592            CALL wad_spg( zsshp2_e, zcpx, zcpy )   ! Calculating W/D gravity filters 
     593            zu_spg(2:jpim1,2:jpjm1) = zu_spg(2:jpim1,2:jpjm1) * zcpx(2:jpim1,2:jpjm1) 
     594            zv_spg(2:jpim1,2:jpjm1) = zv_spg(2:jpim1,2:jpjm1) * zcpy(2:jpim1,2:jpjm1) 
    991595         ENDIF 
    992596         ! 
     
    994598         ! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated 
    995599         ! 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 
     600         ! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated) 
     601         CALL dyn_cor_2d( zhtp2_e, zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV,    zu_trd, zv_trd   ) 
    1057602         ! 
    1058603         ! Add tidal astronomical forcing if defined 
     
    1060605            DO jj = 2, jpjm1 
    1061606               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 
     607                  zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 
     608                  zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 
    1066609               END DO 
    1067610            END DO 
     
    1077620               END DO 
    1078621            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  
     622         ENDIF 
    1104623         ! 
    1105624         ! Set next velocities: 
     625         !     Compute barotropic speeds at step jit+1    (h : total height of the water colomn) 
     626         !--                              VECTOR FORM 
     627         !--     m+1                    m               /                                     m+1/2           \     --! 
     628         !--    u     =                u   + delta_t' * \ grad_x( ssh') -         f * k vect u      +     frc /     --! 
     629         !--                                                                                                        --! 
     630         !--                             FLUX FORM                                                                  --! 
     631         !--   m+1    ___1___   /  m    m               /                  m+1/2              m+1/2    n      \  \  --! 
     632         !--  u    =     m+1   |  h  * u   + delta_t' * \ grad_x( ssh') - h     * f * k vect u      + h * frc /   | --! 
     633         !--            h       \                                                                                /  --! 
     634         !------------------------------------------------------------------------------------------------------------! 
    1106635         IF( ln_dynadv_vec .OR. ln_linssh ) THEN      !* Vector form 
    1107636            DO jj = 2, jpjm1 
    1108637               DO ji = fs_2, fs_jpim1   ! vector opt. 
    1109638                  ua_e(ji,jj) = (                                 un_e(ji,jj)   &  
    1110                             &     + rdtbt * (                      zwx(ji,jj)   & 
     639                            &     + rdtbt * (                   zu_spg(ji,jj)   & 
    1111640                            &                                 + zu_trd(ji,jj)   & 
    1112641                            &                                 + zu_frc(ji,jj) ) &  
     
    1114643 
    1115644                  va_e(ji,jj) = (                                 vn_e(ji,jj)   & 
    1116                             &     + rdtbt * (                      zwy(ji,jj)   & 
     645                            &     + rdtbt * (                   zv_spg(ji,jj)   & 
    1117646                            &                                 + zv_trd(ji,jj)   & 
    1118647                            &                                 + zv_frc(ji,jj) ) & 
    1119648                            &   ) * ssvmask(ji,jj) 
    1120   
    1121649               END DO 
    1122650            END DO 
     
    1125653            DO jj = 2, jpjm1 
    1126654               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 
     655                  !                    ! backward extrapolated depth used in spg terms at jn+1/2 
     656                  zhu_bck = hu_0(ji,jj) + r1_2*r1_e1e2u(ji,jj) * (  e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)    & 
     657                       &                                          + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj)  ) * ssumask(ji,jj) 
     658                  zhv_bck = hv_0(ji,jj) + r1_2*r1_e1e2v(ji,jj) * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )    & 
     659                       &                                          + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1)  ) * ssvmask(ji,jj) 
     660                  !                    ! inverse depth at jn+1 
     661                  z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) ) 
     662                  z1_hv = ssvmask(ji,jj) / ( hu_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) ) 
     663                  ! 
     664                  ua_e(ji,jj) = (               hu_e  (ji,jj) *   un_e (ji,jj)      &  
     665                       &            + rdtbt * (  zhu_bck        * zu_spg (ji,jj)  &   ! 
     666                       &                       + zhup2_e(ji,jj) * zu_trd (ji,jj)  &   ! 
     667                       &                       +  hu_n  (ji,jj) * zu_frc (ji,jj)  )   ) * z1_hu 
     668                  ! 
     669                  va_e(ji,jj) = (               hv_e  (ji,jj) *   vn_e (ji,jj)      & 
     670                       &            + rdtbt * (  zhv_bck        * zv_spg (ji,jj)  &   ! 
     671                       &                       + zhvp2_e(ji,jj) * zv_trd (ji,jj)  &   ! 
     672                       &                       +  hv_n  (ji,jj) * zv_frc (ji,jj)  )   ) * z1_hu 
    1145673               END DO 
    1146674            END DO 
     
    1213741      ! Set advection velocity correction: 
    1214742      IF (ln_bt_fw) THEN 
    1215          zwx(:,:) = un_adv(:,:) 
    1216          zwy(:,:) = vn_adv(:,:) 
    1217743         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(:,:)) 
     744            DO jj = 1, jpj 
     745               DO ji = 1, jpi 
     746                  zun_save = un_adv(ji,jj) 
     747                  zvn_save = vn_adv(ji,jj) 
     748                  !                          ! apply the previously computed correction  
     749                  un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - atfp * un_bf(ji,jj) ) 
     750                  vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - atfp * vn_bf(ji,jj) ) 
     751                  !                          ! Update corrective fluxes for next time step 
     752                  un_bf(ji,jj)  = atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) ) 
     753                  vn_bf(ji,jj)  = atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) ) 
     754                  !                          ! Save integrated transport for next computation 
     755                  ub2_b(ji,jj) = zun_save 
     756                  vb2_b(ji,jj) = zvn_save 
     757               END DO 
     758            END DO 
    1224759         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(:,:) 
     760            un_bf(:,:) = 0._wp            ! corrective fluxes for next time step set to zero 
     761            vn_bf(:,:) = 0._wp 
     762            ub2_b(:,:) = un_adv(:,:)      ! Save integrated transport for next computation 
     763            vb2_b(:,:) = vn_adv(:,:) 
     764         END IF 
    1231765      ENDIF 
    1232766 
     
    14731007      REAL(wp) ::   zxr2, zyr2, zcmax   ! local scalar 
    14741008      REAL(wp), DIMENSION(jpi,jpj) ::   zcu 
     1009      INTEGER  :: inum 
    14751010      !!---------------------------------------------------------------------- 
    14761011      ! 
     
    15791114   END SUBROUTINE dyn_spg_ts_init 
    15801115 
     1116    
     1117   SUBROUTINE dyn_cor_2d_init 
     1118      !!--------------------------------------------------------------------- 
     1119      !!                   ***  ROUTINE dyn_cor_2d_init  *** 
     1120      !! 
     1121      !! ** Purpose : Set time splitting options 
     1122      !! Set arrays to remove/compute coriolis trend. 
     1123      !! Do it once during initialization if volume is fixed, else at each long time step. 
     1124      !! Note that these arrays are also used during barotropic loop. These are however frozen 
     1125      !! although they should be updated in the variable volume case. Not a big approximation. 
     1126      !! To remove this approximation, copy lines below inside barotropic loop 
     1127      !! and update depths at T-F points (ht and zhf resp.) at each barotropic time step 
     1128      !! 
     1129      !! Compute zwz = f / ( height of the water colomn ) 
     1130      !!---------------------------------------------------------------------- 
     1131      INTEGER  ::   ji ,jj, jk              ! dummy loop indices 
     1132      REAL(wp) ::   z1_ht 
     1133      REAL(wp), DIMENSION(jpi,jpj) :: zhf 
     1134      !!---------------------------------------------------------------------- 
     1135      ! 
     1136      ! 
     1137      SELECT CASE( nvor_scheme ) 
     1138      CASE( np_EEN )                != EEN scheme using e3f (energy & enstrophy scheme) 
     1139         SELECT CASE( nn_een_e3f )              !* ff_f/e3 at F-point 
     1140         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4) 
     1141            DO jj = 1, jpjm1 
     1142               DO ji = 1, jpim1 
     1143                  zwz(ji,jj) =   ( ht_n(ji  ,jj+1) + ht_n(ji+1,jj+1) +                    & 
     1144                       &             ht_n(ji  ,jj  ) + ht_n(ji+1,jj  )   ) * 0.25_wp   
     1145                  IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 
     1146               END DO 
     1147            END DO 
     1148         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask) 
     1149            DO jj = 1, jpjm1 
     1150               DO ji = 1, jpim1 
     1151                  zwz(ji,jj) =             (  ht_n  (ji  ,jj+1) + ht_n  (ji+1,jj+1)      & 
     1152                       &                      + ht_n  (ji  ,jj  ) + ht_n  (ji+1,jj  )  )   & 
     1153                       &       / ( MAX( 1._wp,  ssmask(ji  ,jj+1) + ssmask(ji+1,jj+1)      & 
     1154                       &                      + ssmask(ji  ,jj  ) + ssmask(ji+1,jj  )  )   ) 
     1155                  IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 
     1156               END DO 
     1157            END DO 
     1158         END SELECT 
     1159         CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp ) 
     1160         ! 
     1161         ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 
     1162         DO jj = 2, jpj 
     1163            DO ji = 2, jpi 
     1164               ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1) 
     1165               ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  ) 
     1166               ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1) 
     1167               ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  ) 
     1168            END DO 
     1169         END DO 
     1170         ! 
     1171      CASE( np_EET )                  != EEN scheme using e3t (energy conserving scheme) 
     1172         ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 
     1173         DO jj = 2, jpj 
     1174            DO ji = 2, jpi 
     1175               z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) ) 
     1176               ftne(ji,jj) = ( ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) ) * z1_ht 
     1177               ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) ) * z1_ht 
     1178               ftse(ji,jj) = ( ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht 
     1179               ftsw(ji,jj) = ( ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) ) * z1_ht 
     1180            END DO 
     1181         END DO 
     1182         ! 
     1183      CASE( np_ENE, np_ENS , np_MIX )  != all other schemes (ENE, ENS, MIX) except ENT ! 
     1184         ! 
     1185         zwz(:,:) = 0._wp 
     1186         zhf(:,:) = 0._wp 
     1187          
     1188         !!gm  assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed  
     1189!!gm    A priori a better value should be something like : 
     1190!!gm          zhf(i,j) = masked sum of  ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)  
     1191!!gm                     divided by the sum of the corresponding mask  
     1192!!gm  
     1193!!             
     1194         IF( .NOT.ln_sco ) THEN 
     1195   
     1196   !!gm  agree the JC comment  : this should be done in a much clear way 
     1197   
     1198   ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 
     1199   !     Set it to zero for the time being  
     1200   !              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level 
     1201   !              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth 
     1202   !              ENDIF 
     1203   !              zhf(:,:) = gdepw_0(:,:,jk+1) 
     1204            ! 
     1205         ELSE 
     1206            ! 
     1207            !zhf(:,:) = hbatf(:,:) 
     1208            DO jj = 1, jpjm1 
     1209               DO ji = 1, jpim1 
     1210                  zhf(ji,jj) =    (   ht_0  (ji,jj  ) + ht_0  (ji+1,jj  )          & 
     1211                       &            + ht_0  (ji,jj+1) + ht_0  (ji+1,jj+1)   )      & 
     1212                       &     / MAX(   ssmask(ji,jj  ) + ssmask(ji+1,jj  )          & 
     1213                       &            + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp  ) 
     1214               END DO 
     1215            END DO 
     1216         ENDIF 
     1217         ! 
     1218         DO jj = 1, jpjm1 
     1219            zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 
     1220         END DO 
     1221         ! 
     1222         DO jk = 1, jpkm1 
     1223            DO jj = 1, jpjm1 
     1224               zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 
     1225            END DO 
     1226         END DO 
     1227         CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp ) 
     1228         ! JC: TBC. hf should be greater than 0  
     1229         DO jj = 1, jpj 
     1230            DO ji = 1, jpi 
     1231               IF( zhf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zhf(ji,jj) 
     1232            END DO 
     1233         END DO 
     1234         zwz(:,:) = ff_f(:,:) * zwz(:,:) 
     1235      END SELECT 
     1236       
     1237   END SUBROUTINE dyn_cor_2d_init 
     1238 
     1239 
     1240 
     1241   SUBROUTINE dyn_cor_2d( ht_n, hu_n, hv_n, un_b, vn_b, zhU, zhV,    zu_trd, zv_trd   ) 
     1242      !!--------------------------------------------------------------------- 
     1243      !!                   ***  ROUTINE dyn_cor_2d  *** 
     1244      !! 
     1245      !! ** Purpose : Compute u and v coriolis trends 
     1246      !!---------------------------------------------------------------------- 
     1247      INTEGER  ::   ji ,jj               ! dummy loop indices 
     1248      REAL(wp) ::   zx1, zx2, zy1, zy2   !   -      - 
     1249      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: ht_n, hu_n, hv_n, un_b, vn_b, zhU, zhV 
     1250      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: zu_trd, zv_trd 
     1251      !!---------------------------------------------------------------------- 
     1252      SELECT CASE( nvor_scheme ) 
     1253      CASE( np_ENT )                ! enstrophy conserving scheme (f-point) 
     1254         DO jj = 2, jpjm1 
     1255            DO ji = 2, jpim1 
     1256               zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * r1_hu_n(ji,jj)                    & 
     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) * r1_hv_n(ji,jj)                    & 
     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 drg_init( pu_RHSi, pv_RHSi, pCdU_u, pCdU_v ) 
     1451      !!---------------------------------------------------------------------- 
     1452      !!                  ***  ROUTINE 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 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.