Ignore:
Timestamp:
2019-11-22T15:29:17+01:00 (20 months ago)
Author:
acc
Message:

Merge in changes from 2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps. This just creates a fresh copy of this branch to use as the merge base. See ticket #2341

Location:
NEMO/branches/2019/dev_r11943_MERGE_2019/src
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11943_MERGE_2019/src

    • Property svn:mergeinfo deleted
  • NEMO/branches/2019/dev_r11943_MERGE_2019/src/NST/agrif_oce_update.F90

    r10068 r11949  
    230230      ! ----------------------- 
    231231      ! 
    232       e3u_a(:,:,:) = e3u_n(:,:,:) 
    233       e3v_a(:,:,:) = e3v_n(:,:,:) 
    234 !      ua(:,:,:) = e3u_b(:,:,:) 
    235 !      va(:,:,:) = e3v_b(:,:,:) 
    236       hu_a(:,:) = hu_n(:,:) 
    237       hv_a(:,:) = hv_n(:,:) 
     232      e3u(:,:,:,Krhs_a) = e3u(:,:,:,Kmm_a) 
     233      e3v(:,:,:,Krhs_a) = e3v(:,:,:,Kmm_a) 
     234!      uu(:,:,:,Krhs_a) = e3u(:,:,:,Kbb_a) 
     235!      vv(:,:,:,Krhs_a) = e3v(:,:,:,Kbb_a) 
     236      hu(:,:,Krhs_a) = hu(:,:,Kmm_a) 
     237      hv(:,:,Krhs_a) = hv(:,:,Kmm_a) 
    238238 
    239239      ! 1) NOW fields 
     
    242242         ! Vertical scale factor interpolations 
    243243         ! ------------------------------------ 
    244       CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:) ,  'U' ) 
    245       CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:) ,  'V' ) 
    246       CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:) ,  'F' ) 
    247  
    248       CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 
    249       CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 
     244      CALL dom_vvl_interpol( e3t(:,:,:,Kmm_a), e3u(:,:,:,Kmm_a) ,  'U' ) 
     245      CALL dom_vvl_interpol( e3t(:,:,:,Kmm_a), e3v(:,:,:,Kmm_a) ,  'V' ) 
     246      CALL dom_vvl_interpol( e3u(:,:,:,Kmm_a), e3f(:,:,:) ,  'F' ) 
     247 
     248      CALL dom_vvl_interpol( e3u(:,:,:,Kmm_a), e3uw(:,:,:,Kmm_a), 'UW' ) 
     249      CALL dom_vvl_interpol( e3v(:,:,:,Kmm_a), e3vw(:,:,:,Kmm_a), 'VW' ) 
    250250 
    251251         ! Update total depths: 
    252252         ! -------------------- 
    253       hu_n(:,:) = 0._wp                        ! Ocean depth at U-points 
    254       hv_n(:,:) = 0._wp                        ! Ocean depth at V-points 
     253      hu(:,:,Kmm_a) = 0._wp                        ! Ocean depth at U-points 
     254      hv(:,:,Kmm_a) = 0._wp                        ! Ocean depth at V-points 
    255255      DO jk = 1, jpkm1 
    256          hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk) 
    257          hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk) 
     256         hu(:,:,Kmm_a) = hu(:,:,Kmm_a) + e3u(:,:,jk,Kmm_a) * umask(:,:,jk) 
     257         hv(:,:,Kmm_a) = hv(:,:,Kmm_a) + e3v(:,:,jk,Kmm_a) * vmask(:,:,jk) 
    258258      END DO 
    259259      !                                        ! Inverse of the local depth 
    260       r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) 
    261       r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) 
     260      r1_hu(:,:,Kmm_a) = ssumask(:,:) / ( hu(:,:,Kmm_a) + 1._wp - ssumask(:,:) ) 
     261      r1_hv(:,:,Kmm_a) = ssvmask(:,:) / ( hv(:,:,Kmm_a) + 1._wp - ssvmask(:,:) ) 
    262262 
    263263 
     
    268268         ! Vertical scale factor interpolations 
    269269         ! ------------------------------------ 
    270          CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:),  'U'  ) 
    271          CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:),  'V'  ) 
    272  
    273          CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 
    274          CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 
     270         CALL dom_vvl_interpol( e3t(:,:,:,Kbb_a), e3u(:,:,:,Kbb_a),  'U'  ) 
     271         CALL dom_vvl_interpol( e3t(:,:,:,Kbb_a), e3v(:,:,:,Kbb_a),  'V'  ) 
     272 
     273         CALL dom_vvl_interpol( e3u(:,:,:,Kbb_a), e3uw(:,:,:,Kbb_a), 'UW' ) 
     274         CALL dom_vvl_interpol( e3v(:,:,:,Kbb_a), e3vw(:,:,:,Kbb_a), 'VW' ) 
    275275 
    276276         ! Update total depths: 
    277277         ! -------------------- 
    278          hu_b(:,:) = 0._wp                     ! Ocean depth at U-points 
    279          hv_b(:,:) = 0._wp                     ! Ocean depth at V-points 
     278         hu(:,:,Kbb_a) = 0._wp                     ! Ocean depth at U-points 
     279         hv(:,:,Kbb_a) = 0._wp                     ! Ocean depth at V-points 
    280280         DO jk = 1, jpkm1 
    281             hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk) 
    282             hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk) 
     281            hu(:,:,Kbb_a) = hu(:,:,Kbb_a) + e3u(:,:,jk,Kbb_a) * umask(:,:,jk) 
     282            hv(:,:,Kbb_a) = hv(:,:,Kbb_a) + e3v(:,:,jk,Kbb_a) * vmask(:,:,jk) 
    283283         END DO 
    284284         !                                     ! Inverse of the local depth 
    285          r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) 
    286          r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) 
     285         r1_hu(:,:,Kbb_a) = ssumask(:,:) / ( hu(:,:,Kbb_a) + 1._wp - ssumask(:,:) ) 
     286         r1_hv(:,:,Kbb_a) = ssvmask(:,:) / ( hv(:,:,Kbb_a) + 1._wp - ssvmask(:,:) ) 
    287287      ENDIF 
    288288      ! 
     
    315315               DO jj=j1,j2 
    316316                  DO ji=i1,i2 
    317                      tabres(ji,jj,jk,jn) = (tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) ) & 
    318                                            * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp 
     317                     tabres(ji,jj,jk,jn) = (    ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) )          & 
     318                       &                   * tmask(ji,jj,jk)       + (tmask(ji,jj,jk) - 1._wp)*999._wp 
    319319                  END DO 
    320320               END DO 
     
    324324            DO jj=j1,j2 
    325325               DO ji=i1,i2 
    326                   tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) & 
    327                                            + (tmask(ji,jj,jk)-1)*999._wp 
     326                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) & 
     327                    &                      + (tmask(ji,jj,jk) - 1._wp)*999._wp 
    328328               END DO 
    329329            END DO 
    330330         END DO 
    331331      ELSE 
    332          tabres_child(:,:,:,:) = 0. 
     332         tabres_child(:,:,:,:) = 0._wp 
    333333         AGRIF_SpecialValue = 0._wp 
    334334         DO jj=j1,j2 
     
    345345                  IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF 
    346346                  N_out = N_out + 1 
    347                   h_out(N_out) = e3t_n(ji,jj,jk)  
     347                  h_out(N_out) = e3t(ji,jj,jk,Kmm_a)  
    348348               ENDDO 
    349349               IF (N_in > 0) THEN !Remove this? 
     
    356356                  ENDIF 
    357357                  DO jn=n1,n2-1 
    358                      CALL reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),tabres_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out) 
     358                     CALL reconstructandremap( tabin(1:N_in,jn), h_in(1:N_in), tabres_child(ji,jj,1:N_out,jn),  & 
     359                       &                       h_out(1:N_out)  , N_in        , N_out                           ) 
    359360                  ENDDO 
    360361               ENDIF 
     
    369370                     DO ji=i1,i2 
    370371                        IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 
    371                            tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) &  
    372                                  & + atfp * ( tabres_child(ji,jj,jk,jn) & 
    373                                  &          - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 
     372                           ts(ji,jj,jk,jn,Kbb_a) =                     ts(ji,jj,jk,jn,Kbb_a) &  
     373                                 &                + atfp * ( tabres_child(ji,jj,jk,jn)       & 
     374                                 &                                   - ts(ji,jj,jk,jn,Kmm_a) & 
     375                                 &                         ) * tmask(ji,jj,jk) 
    374376                        ENDIF 
    375377                     ENDDO 
     
    383385                  DO ji=i1,i2 
    384386                     IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN  
    385                         tsn(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk) 
     387                        ts(ji,jj,jk,jn,Kmm_a) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk) 
    386388                     END IF 
    387389                  END DO 
     
    413415                  DO ji=i1,i2 
    414416!> jc tmp 
    415                      tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk) / e3t_0(ji,jj,jk) 
    416 !                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk) 
     417                     tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a)  * e3t(ji,jj,jk,Kmm_a) / e3t_0(ji,jj,jk) 
     418!                    tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a)  * e3t(ji,jj,jk,Kmm_a) 
    417419!< jc tmp 
    418420                  END DO 
     
    434436                     DO ji = i1, i2 
    435437                        IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN 
    436                            ztb  = tsb(ji,jj,jk,jn) * e3t_b(ji,jj,jk) ! fse3t_b prior update should be used 
     438                           ztb  = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 
    437439                           ztnu = tabres(ji,jj,jk,jn) 
    438                            ztno = tsn(ji,jj,jk,jn) * e3t_a(ji,jj,jk) 
    439                            tsb(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  &  
    440                                      &        * tmask(ji,jj,jk) / e3t_b(ji,jj,jk) 
     440                           ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a) 
     441                           ts(ji,jj,jk,jn,Kbb_a) = ( ztb + atfp * ( ztnu - ztno) )  &  
     442                                     &            * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a) 
    441443                        ENDIF 
    442444                     END DO 
     
    450452                  DO ji=i1,i2 
    451453                     IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN  
    452                         tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) / e3t_n(ji,jj,jk) 
     454                        ts(ji,jj,jk,jn,Kmm_a) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm_a) 
    453455                     END IF 
    454456                  END DO 
     
    458460         ! 
    459461         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
    460             tsb(i1:i2,j1:j2,k1:k2,1:jpts)  = tsn(i1:i2,j1:j2,k1:k2,1:jpts) 
     462            ts(i1:i2,j1:j2,k1:k2,1:jpts,Kbb_a)  = ts(i1:i2,j1:j2,k1:k2,1:jpts,Kmm_a) 
    461463         ENDIF 
    462464         ! 
     
    495497            DO jj=j1,j2 
    496498               DO ji=i1,i2 
    497                   tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk)  & 
    498                                        + (umask(ji,jj,jk)-1)*999._wp 
    499                   tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)  & 
    500                                        + (umask(ji,jj,jk)-1)*999._wp 
     499                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * umask(ji,jj,jk) * uu(ji,jj,jk,Kmm_a)  & 
     500                    &                  + (umask(ji,jj,jk)-1._wp)*999._wp 
     501                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)  & 
     502                    &                  + (umask(ji,jj,jk)-1._wp)*999._wp 
    501503               END DO 
    502504            END DO 
     
    513515                  IF( tabres(ji,jj,jk,2) < -900) EXIT 
    514516                  N_in = N_in + 1 
    515                   tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
     517                  tabin(jk)  = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
    516518                  h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj) 
    517519               ENDDO 
     
    520522                  IF (umask(ji,jj,jk) == 0) EXIT 
    521523                  N_out = N_out + 1 
    522                   h_out(N_out) = e3u_n(ji,jj,jk) 
     524                  h_out(N_out) = e3u(ji,jj,jk,Kmm_a) 
    523525               ENDDO 
    524526               IF (N_in * N_out > 0) THEN 
    525                   h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     527                  h_diff = SUM( h_out(1:N_out) ) - SUM( h_in(1:N_in) ) 
    526528                  IF (h_diff < -1.e-4) THEN 
    527529!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.  
     
    540542                     ENDDO 
    541543                  ENDIF 
    542                   CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 
    543                   tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out)) 
     544                  CALL reconstructandremap( tabin(1:N_in) , h_in(1:N_in), tabres_child(ji,jj,1:N_out),       & 
     545                    &                       h_out(1:N_out), N_in        , N_out                            ) 
     546                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/( e2u(ji,jj)*h_out(N_out) ) 
    544547               ENDIF 
    545548            ENDDO 
     
    550553               DO ji=i1,i2 
    551554                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    552                      ub(ji,jj,jk) = ub(ji,jj,jk) &  
    553                            & + atfp * ( tabres_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk) 
     555                     uu(ji,jj,jk,Kbb_a) =                     uu(ji,jj,jk,Kbb_a)                                  &  
     556                           &             + atfp * ( tabres_child(ji,jj,jk) - uu(ji,jj,jk,Kmm_a) ) * umask(ji,jj,jk) 
    554557                  ENDIF 
    555558                  ! 
    556                   un(ji,jj,jk) = tabres_child(ji,jj,jk) * umask(ji,jj,jk) 
     559                  uu(ji,jj,jk,Kmm_a) = tabres_child(ji,jj,jk) * umask(ji,jj,jk) 
    557560               END DO 
    558561            END DO 
     
    579582         zrhoy = Agrif_Rhoy() 
    580583         DO jk = k1, k2 
    581             tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk) 
     584            tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm_a) & 
     585              &                                                 *  uu(i1:i2,j1:j2,jk,Kmm_a) 
    582586         END DO 
    583587      ELSE 
     
    588592                  ! 
    589593                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    590                      zub  = ub(ji,jj,jk) * e3u_b(ji,jj,jk)  ! fse3t_b prior update should be used 
    591                      zuno = un(ji,jj,jk) * e3u_a(ji,jj,jk) 
     594                     zub  = uu(ji,jj,jk,Kbb_a) * e3u(ji,jj,jk,Kbb_a)  ! fse3t_b prior update should be used 
     595                     zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a) 
    592596                     zunu = tabres(ji,jj,jk,1) 
    593                      ub(ji,jj,jk) = ( zub + atfp * ( zunu - zuno) ) &       
    594                                     & * umask(ji,jj,jk) / e3u_b(ji,jj,jk) 
     597                     uu(ji,jj,jk,Kbb_a) = ( zub + atfp * ( zunu - zuno ) )      &       
     598                       &                  * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a) 
    595599                  ENDIF 
    596600                  ! 
    597                   un(ji,jj,jk) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u_n(ji,jj,jk) 
     601                  uu(ji,jj,jk,Kmm_a) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u(ji,jj,jk,Kmm_a) 
    598602               END DO 
    599603            END DO 
     
    601605         ! 
    602606         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
    603             ub(i1:i2,j1:j2,k1:k2)  = un(i1:i2,j1:j2,k1:k2) 
     607            uu(i1:i2,j1:j2,k1:k2,Kbb_a)  = uu(i1:i2,j1:j2,k1:k2,Kmm_a) 
    604608         ENDIF 
    605609         ! 
     
    632636         IF (western_side) THEN 
    633637            DO jj=j1,j2 
    634                zcor = un_b(i1-1,jj) * hu_a(i1-1,jj) * r1_hu_n(i1-1,jj) - un_b(i1-1,jj) 
    635                un_b(i1-1,jj) = un_b(i1-1,jj) + zcor 
     638               zcor = uu_b(i1-1,jj,Kmm_a) * hu(i1-1,jj,Krhs_a) * r1_hu(i1-1,jj,Kmm_a) - uu_b(i1-1,jj,Kmm_a) 
     639               uu_b(i1-1,jj,Kmm_a) = uu_b(i1-1,jj,Kmm_a) + zcor 
    636640               DO jk=1,jpkm1 
    637                   un(i1-1,jj,jk) = un(i1-1,jj,jk) + zcor * umask(i1-1,jj,jk) 
     641                  uu(i1-1,jj,jk,Kmm_a) = uu(i1-1,jj,jk,Kmm_a) + zcor * umask(i1-1,jj,jk) 
    638642               END DO  
    639643            END DO 
     
    642646         IF (eastern_side) THEN 
    643647            DO jj=j1,j2 
    644                zcor = un_b(i2+1,jj) * hu_a(i2+1,jj) * r1_hu_n(i2+1,jj) - un_b(i2+1,jj) 
    645                un_b(i2+1,jj) = un_b(i2+1,jj) + zcor 
     648               zcor = uu_b(i2+1,jj,Kmm_a) * hu(i2+1,jj,Krhs_a) * r1_hu(i2+1,jj,Kmm_a) - uu_b(i2+1,jj,Kmm_a) 
     649               uu_b(i2+1,jj,Kmm_a) = uu_b(i2+1,jj,Kmm_a) + zcor 
    646650               DO jk=1,jpkm1 
    647                   un(i2+1,jj,jk) = un(i2+1,jj,jk) + zcor * umask(i2+1,jj,jk) 
     651                  uu(i2+1,jj,jk,Kmm_a) = uu(i2+1,jj,jk,Kmm_a) + zcor * umask(i2+1,jj,jk) 
    648652               END DO  
    649653            END DO 
     
    682686            DO jj=j1,j2 
    683687               DO ji=i1,i2 
    684                   tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) & 
    685                                        + (vmask(ji,jj,jk)-1)*999._wp 
    686                   tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) & 
     688                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a)                   & 
     689                    &                  *  vmask(ji,jj,jk) *  vv(ji,jj,jk,Kmm_a)                   & 
     690                    &                  + (vmask(ji,jj,jk)-1)*999._wp 
     691                  tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) & 
    687692                                       + (vmask(ji,jj,jk)-1)*999._wp 
    688693               END DO 
     
    705710                  IF (vmask(ji,jj,jk) == 0) EXIT 
    706711                  N_out = N_out + 1 
    707                   h_out(N_out) = e3v_n(ji,jj,jk) 
     712                  h_out(N_out) = e3v(ji,jj,jk,Kmm_a) 
    708713               ENDDO 
    709714               IF (N_in * N_out > 0) THEN 
    710                   h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     715                  h_diff = SUM( h_out(1:N_out) ) - SUM( h_in(1:N_in) ) 
    711716                  IF (h_diff < -1.e-4) then 
    712717!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.  
     
    714719                     excess = 0._wp 
    715720                     DO jk=N_in,1,-1 
    716                         thick = MIN(-1*h_diff, h_in(jk)) 
     721                        thick = MIN( -1._wp*h_diff, h_in(jk) ) 
    717722                        excess = excess + tabin(jk)*thick*e2u(ji,jj) 
    718                         tabin(jk) = tabin(jk)*(1. - thick/h_in(jk)) 
     723                        tabin(jk) = tabin(jk)*(1._wp - thick/h_in(jk)) 
    719724                        h_diff = h_diff + thick 
    720725                        IF ( h_diff == 0) THEN 
     
    725730                     ENDDO 
    726731                  ENDIF 
    727                   CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 
     732                  CALL reconstructandremap( tabin(1:N_in) , h_in(1:N_in), tabres_child(ji,jj,1:N_out),     & 
     733                    &                       h_out(1:N_out), N_in        , N_out                          ) 
    728734                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out)) 
    729735               ENDIF 
     
    736742                  ! 
    737743                  IF( .NOT.(lk_agrif_fstep.AND.(neuler==0)) ) THEN ! Add asselin part 
    738                      vb(ji,jj,jk) = vb(ji,jj,jk) &  
    739                            & + atfp * ( tabres_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk) 
     744                     vv(ji,jj,jk,Kbb_a) =                     vv(ji,jj,jk,Kbb_a)                  &  
     745                           &             + atfp * ( tabres_child(ji,jj,jk) - vv(ji,jj,jk,Kmm_a) ) & 
     746                           &                      *        vmask(ji,jj,jk) 
    740747                  ENDIF 
    741748                  ! 
    742                   vn(ji,jj,jk) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk) 
     749                  vv(ji,jj,jk,Kmm_a) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk) 
    743750               END DO 
    744751            END DO 
     
    767774            DO jj=j1,j2 
    768775               DO ji=i1,i2 
    769                   tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk) 
     776                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a) 
    770777               END DO 
    771778            END DO 
     
    778785                  ! 
    779786                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    780                      zvb  = vb(ji,jj,jk) * e3v_b(ji,jj,jk) ! fse3t_b prior update should be used 
    781                      zvno = vn(ji,jj,jk) * e3v_a(ji,jj,jk) 
     787                     zvb  = vv(ji,jj,jk,Kbb_a) * e3v(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 
     788                     zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a) 
    782789                     zvnu = tabres(ji,jj,jk,1) 
    783                      vb(ji,jj,jk) = ( zvb + atfp * ( zvnu - zvno) ) &       
    784                                     & * vmask(ji,jj,jk) / e3v_b(ji,jj,jk) 
     790                     vv(ji,jj,jk,Kbb_a) = ( zvb + atfp * ( zvnu - zvno) )        &       
     791                       &                  * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a) 
    785792                  ENDIF 
    786793                  ! 
    787                   vn(ji,jj,jk) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v_n(ji,jj,jk) 
     794                  vv(ji,jj,jk,Kmm_a) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kmm_a) 
    788795               END DO 
    789796            END DO 
     
    791798         ! 
    792799         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
    793             vb(i1:i2,j1:j2,k1:k2)  = vn(i1:i2,j1:j2,k1:k2) 
     800            vv(i1:i2,j1:j2,k1:k2,Kbb_a)  = vv(i1:i2,j1:j2,k1:k2,Kmm_a) 
    794801         ENDIF 
    795802         ! 
     
    822829         IF (southern_side) THEN 
    823830            DO ji=i1,i2 
    824                zcor = vn_b(ji,j1-1) * hv_a(ji,j1-1) * r1_hv_n(ji,j1-1) - vn_b(ji,j1-1) 
    825                vn_b(ji,j1-1) = vn_b(ji,j1-1) + zcor 
     831               zcor = vv_b(ji,j1-1,Kmm_a) * hv(ji,j1-1,Krhs_a) * r1_hv(ji,j1-1,Kmm_a) - vv_b(ji,j1-1,Kmm_a) 
     832               vv_b(ji,j1-1,Kmm_a) = vv_b(ji,j1-1,Kmm_a) + zcor 
    826833               DO jk=1,jpkm1 
    827                   vn(ji,j1-1,jk) = vn(ji,j1-1,jk) + zcor * vmask(ji,j1-1,jk) 
     834                  vv(ji,j1-1,jk,Kmm_a) = vv(ji,j1-1,jk,Kmm_a) + zcor * vmask(ji,j1-1,jk) 
    828835               END DO  
    829836            END DO 
     
    832839         IF (northern_side) THEN 
    833840            DO ji=i1,i2 
    834                zcor = vn_b(ji,j2+1) * hv_a(ji,j2+1) * r1_hv_n(ji,j2+1) - vn_b(ji,j2+1) 
    835                vn_b(ji,j2+1) = vn_b(ji,j2+1) + zcor 
     841               zcor = vv_b(ji,j2+1,Kmm_a) * hv(ji,j2+1,Krhs_a) * r1_hv(ji,j2+1,Kmm_a) - vv_b(ji,j2+1,Kmm_a) 
     842               vv_b(ji,j2+1,Kmm_a) = vv_b(ji,j2+1,Kmm_a) + zcor 
    836843               DO jk=1,jpkm1 
    837                   vn(ji,j2+1,jk) = vn(ji,j2+1,jk) + zcor * vmask(ji,j2+1,jk) 
     844                  vv(ji,j2+1,jk,Kmm_a) = vv(ji,j2+1,jk,Kmm_a) + zcor * vmask(ji,j2+1,jk) 
    838845               END DO  
    839846            END DO 
     
    862869         DO jj=j1,j2 
    863870            DO ji=i1,i2 
    864                tabres(ji,jj) = zrhoy * un_b(ji,jj) * hu_n(ji,jj) * e2u(ji,jj) 
     871               tabres(ji,jj) = zrhoy * uu_b(ji,jj,Kmm_a) * hu(ji,jj,Kmm_a) * e2u(ji,jj) 
    865872            END DO 
    866873         END DO 
     
    873880               spgu(ji,jj) = 0._wp 
    874881               DO jk=1,jpkm1 
    875                   spgu(ji,jj) = spgu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) 
     882                  spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a) 
    876883               END DO 
    877884               ! 
    878                zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu_n(ji,jj) 
     885               zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu(ji,jj,Kmm_a) 
    879886               DO jk=1,jpkm1               
    880                   un(ji,jj,jk) = un(ji,jj,jk) + zcorr * umask(ji,jj,jk)            
     887                  uu(ji,jj,jk,Kmm_a) = uu(ji,jj,jk,Kmm_a) + zcorr * umask(ji,jj,jk)            
    881888               END DO 
    882889               ! 
    883890               ! Update barotropic velocities: 
    884891               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 
    885                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    886                      zcorr = (tabres(ji,jj) - un_b(ji,jj) * hu_a(ji,jj)) * r1_hu_b(ji,jj) 
    887                      ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1) 
     892                  IF ( .NOT.( lk_agrif_fstep .AND. (neuler==0) ) ) THEN                          ! Add asselin part 
     893                     zcorr = (tabres(ji,jj) - uu_b(ji,jj,Kmm_a) * hu(ji,jj,Krhs_a)) * r1_hu(ji,jj,Kbb_a) 
     894                     uu_b(ji,jj,Kbb_a) = uu_b(ji,jj,Kbb_a) + atfp * zcorr * umask(ji,jj,1) 
    888895                  END IF 
    889896               ENDIF     
    890                un_b(ji,jj) = tabres(ji,jj) * r1_hu_n(ji,jj) * umask(ji,jj,1) 
     897               uu_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hu(ji,jj,Kmm_a) * umask(ji,jj,1) 
    891898               !        
    892899               ! Correct "before" velocities to hold correct bt component: 
    893900               spgu(ji,jj) = 0.e0 
    894901               DO jk=1,jpkm1 
    895                   spgu(ji,jj) = spgu(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) 
     902                  spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kbb_a) * uu(ji,jj,jk,Kbb_a) 
    896903               END DO 
    897904               ! 
    898                zcorr = ub_b(ji,jj) - spgu(ji,jj) * r1_hu_b(ji,jj) 
     905               zcorr = uu_b(ji,jj,Kbb_a) - spgu(ji,jj) * r1_hu(ji,jj,Kbb_a) 
    899906               DO jk=1,jpkm1               
    900                   ub(ji,jj,jk) = ub(ji,jj,jk) + zcorr * umask(ji,jj,jk)            
     907                  uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) + zcorr * umask(ji,jj,jk)            
    901908               END DO 
    902909               ! 
     
    905912         ! 
    906913         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
    907             ub_b(i1:i2,j1:j2)  = un_b(i1:i2,j1:j2) 
     914            uu_b(i1:i2,j1:j2,Kbb_a)  = uu_b(i1:i2,j1:j2,Kmm_a) 
    908915         ENDIF 
    909916      ENDIF 
     
    928935         DO jj=j1,j2 
    929936            DO ji=i1,i2 
    930                tabres(ji,jj) = zrhox * vn_b(ji,jj) * hv_n(ji,jj) * e1v(ji,jj)  
     937               tabres(ji,jj) = zrhox * vv_b(ji,jj,Kmm_a) * hv(ji,jj,Kmm_a) * e1v(ji,jj)  
    931938            END DO 
    932939         END DO 
     
    939946               spgv(ji,jj) = 0.e0 
    940947               DO jk=1,jpkm1 
    941                   spgv(ji,jj) = spgv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) 
     948                  spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a) 
    942949               END DO 
    943950               ! 
    944                zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv_n(ji,jj) 
     951               zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv(ji,jj,Kmm_a) 
    945952               DO jk=1,jpkm1               
    946                   vn(ji,jj,jk) = vn(ji,jj,jk) + zcorr * vmask(ji,jj,jk)            
     953                  vv(ji,jj,jk,Kmm_a) = vv(ji,jj,jk,Kmm_a) + zcorr * vmask(ji,jj,jk)            
    947954               END DO 
    948955               ! 
    949956               ! Update barotropic velocities: 
    950                IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 
    951                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    952                      zcorr = (tabres(ji,jj) - vn_b(ji,jj) * hv_a(ji,jj)) * r1_hv_b(ji,jj) 
    953                      vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1) 
     957               IF ( .NOT.ln_dynspg_ts .OR. ( ln_dynspg_ts .AND. ( .NOT.ln_bt_fw ) ) ) THEN 
     958                  IF ( .NOT. ( lk_agrif_fstep .AND. ( neuler==0 ) ) ) THEN                      ! Add asselin part 
     959                     zcorr = (tabres(ji,jj) - vv_b(ji,jj,Kmm_a) * hv(ji,jj,Krhs_a)) * r1_hv(ji,jj,Kbb_a) 
     960                     vv_b(ji,jj,Kbb_a) =       vv_b(ji,jj,Kbb_a) + atfp * zcorr  *  vmask(ji,jj,1) 
    954961                  END IF 
    955962               ENDIF               
    956                vn_b(ji,jj) = tabres(ji,jj) * r1_hv_n(ji,jj) * vmask(ji,jj,1) 
     963               vv_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hv(ji,jj,Kmm_a) * vmask(ji,jj,1) 
    957964               !        
    958965               ! Correct "before" velocities to hold correct bt component: 
    959966               spgv(ji,jj) = 0.e0 
    960967               DO jk=1,jpkm1 
    961                   spgv(ji,jj) = spgv(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) 
     968                  spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kbb_a) * vv(ji,jj,jk,Kbb_a) 
    962969               END DO 
    963970               ! 
    964                zcorr = vb_b(ji,jj) - spgv(ji,jj) * r1_hv_b(ji,jj) 
     971               zcorr = vv_b(ji,jj,Kbb_a) - spgv(ji,jj) * r1_hv(ji,jj,Kbb_a) 
    965972               DO jk=1,jpkm1               
    966                   vb(ji,jj,jk) = vb(ji,jj,jk) + zcorr * vmask(ji,jj,jk)            
     973                  vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) + zcorr * vmask(ji,jj,jk)            
    967974               END DO 
    968975               ! 
     
    971978         ! 
    972979         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
    973             vb_b(i1:i2,j1:j2)  = vn_b(i1:i2,j1:j2) 
     980            vv_b(i1:i2,j1:j2,Kbb_a)  = vv_b(i1:i2,j1:j2,Kmm_a) 
    974981         ENDIF 
    975982         ! 
     
    9931000         DO jj=j1,j2 
    9941001            DO ji=i1,i2 
    995                tabres(ji,jj) = sshn(ji,jj) 
     1002               tabres(ji,jj) = ssh(ji,jj,Kmm_a) 
    9961003            END DO 
    9971004         END DO 
     
    10001007            DO jj=j1,j2 
    10011008               DO ji=i1,i2 
    1002                   sshb(ji,jj) =   sshb(ji,jj) & 
    1003                         & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1) 
     1009                  ssh(ji,jj,Kbb_a) =              ssh(ji,jj,Kbb_a)                & 
     1010                    &               + atfp * ( tabres(ji,jj) - ssh(ji,jj,Kmm_a) ) & 
     1011                    &                        *  tmask(ji,jj,1) 
    10041012               END DO 
    10051013            END DO 
     
    10081016         DO jj=j1,j2 
    10091017            DO ji=i1,i2 
    1010                sshn(ji,jj) = tabres(ji,jj) * tmask(ji,jj,1) 
     1018               ssh(ji,jj,Kmm_a) = tabres(ji,jj) * tmask(ji,jj,1) 
    10111019            END DO 
    10121020         END DO 
    10131021         ! 
    10141022         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
    1015             sshb(i1:i2,j1:j2)  = sshn(i1:i2,j1:j2) 
     1023            ssh(i1:i2,j1:j2,Kbb_a)  = ssh(i1:i2,j1:j2,Kmm_a) 
    10161024         ENDIF 
    10171025         ! 
     
    10941102            DO jj=j1,j2 
    10951103               zcor = rdt * r1_e1e2t(i1  ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj))  
    1096                sshn(i1  ,jj) = sshn(i1  ,jj) + zcor 
    1097                IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i1  ,jj) = sshb(i1  ,jj) + atfp * zcor 
     1104               ssh(i1  ,jj,Kmm_a) = ssh(i1  ,jj,Kmm_a) + zcor 
     1105               IF ( .NOT. ( lk_agrif_fstep .AND. ( neuler==0 ) ) )  & 
     1106                 &                  ssh(i1  ,jj,Kbb_a) = ssh(i1  ,jj,Kbb_a) + atfp * zcor 
    10981107            END DO 
    10991108         ENDIF 
     
    11011110            DO jj=j1,j2 
    11021111               zcor = - rdt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj)) 
    1103                sshn(i2+1,jj) = sshn(i2+1,jj) + zcor 
    1104                IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i2+1,jj) = sshb(i2+1,jj) + atfp * zcor 
     1112               ssh(i2+1,jj,Kmm_a) = ssh(i2+1,jj,Kmm_a) + zcor 
     1113               IF ( .NOT. ( lk_agrif_fstep .AND. ( neuler==0 ) ) )   & 
     1114                 &                  ssh(i2+1,jj,Kbb_a) = ssh(i2+1,jj,Kbb_a) + atfp * zcor 
    11051115            END DO 
    11061116         ENDIF 
     
    11811191         IF (southern_side) THEN 
    11821192            DO ji=i1,i2 
    1183                zcor = rdt * r1_e1e2t(ji,j1  ) * e1v(ji,j1  ) * (vb2_b(ji,j1)-tabres(ji,j1)) 
    1184                sshn(ji,j1  ) = sshn(ji,j1  ) + zcor 
    1185                IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j1  ) = sshb(ji,j1) + atfp * zcor 
     1193               zcor = rdt * r1_e1e2t(ji,j1  ) * e1v(ji,j1  ) * ( vb2_b(ji,j1)-tabres(ji,j1) ) 
     1194               ssh(ji,j1  ,Kmm_a) = ssh(ji,j1  ,Kmm_a) + zcor 
     1195               IF ( .NOT. ( lk_agrif_fstep .AND. ( neuler==0 ) ) )   & 
     1196                 &                 ssh(ji,j1  ,Kbb_a) = ssh(ji,j1,Kbb_a) + atfp * zcor 
    11861197            END DO 
    11871198         ENDIF 
    11881199         IF (northern_side) THEN                
    11891200            DO ji=i1,i2 
    1190                zcor = - rdt * r1_e1e2t(ji,j2+1) * e1v(ji,j2  ) * (vb2_b(ji,j2)-tabres(ji,j2)) 
    1191                sshn(ji,j2+1) = sshn(ji,j2+1) + zcor 
    1192                IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j2+1) = sshb(ji,j2+1) + atfp * zcor 
     1201               zcor = - rdt * r1_e1e2t(ji,j2+1) * e1v(ji,j2  ) * ( vb2_b(ji,j2)-tabres(ji,j2) ) 
     1202               ssh(ji,j2+1,Kmm_a) = ssh(ji,j2+1,Kmm_a) + zcor 
     1203               IF ( .NOT. ( lk_agrif_fstep .AND. ( neuler==0 ) ) )   & 
     1204                 &                  ssh(ji,j2+1,Kbb_a) = ssh(ji,j2+1,Kbb_a) + atfp * zcor 
    11931205            END DO 
    11941206         ENDIF 
     
    12231235            END DO 
    12241236         END DO 
    1225          tabres(:,:,:,1)=tabres(:,:,:,1)*Agrif_Rhox()*Agrif_Rhoy() 
    1226          tabres(:,:,:,2)=tabres(:,:,:,2)*Agrif_Rhox() 
    1227          tabres(:,:,:,3)=tabres(:,:,:,3)*Agrif_Rhoy() 
     1237         tabres(:,:,:,1) = tabres(:,:,:,1)*Agrif_Rhox()*Agrif_Rhoy() 
     1238         tabres(:,:,:,2) = tabres(:,:,:,2)*Agrif_Rhox() 
     1239         tabres(:,:,:,3) = tabres(:,:,:,3)*Agrif_Rhoy() 
    12281240      ELSE 
    12291241         DO jk=k1,k2 
     
    12341246                     print *,'VAL2 = ',ji,jj,jk,tabres(ji,jj,jk,2),e1t(ji,jj)*tmask(ji,jj,jk) 
    12351247                     print *,'VAL3 = ',ji,jj,jk,tabres(ji,jj,jk,3),e2t(ji,jj)*tmask(ji,jj,jk) 
    1236                      ztemp = sqrt(tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3))) 
     1248                     ztemp = SQRT( tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3)) ) 
    12371249                     print *,'CORR = ',ztemp-1. 
    12381250                     print *,'NEW VALS = ',tabres(ji,jj,jk,2)*ztemp,tabres(ji,jj,jk,3)*ztemp, & 
    12391251                           tabres(ji,jj,jk,2)*ztemp*tabres(ji,jj,jk,3)*ztemp 
    1240                      e1t(ji,jj) = tabres(ji,jj,jk,2)*ztemp 
    1241                      e2t(ji,jj) = tabres(ji,jj,jk,3)*ztemp 
     1252                     e1t(ji,jj) = tabres(ji,jj,jk,2) * ztemp 
     1253                     e2t(ji,jj) = tabres(ji,jj,jk,3) * ztemp 
    12421254                  END IF 
    12431255               END DO 
     
    13191331            DO jj=j1,j2 
    13201332               DO ji=i1,i2 
    1321                   ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + sshn(ji,jj) & 
    1322                                      & *ssmask(ji,jj)/(ht_0(ji,jj)-1._wp + ssmask(ji,jj))) 
     1333                  ptab(ji,jj,jk) =   e3t_0(ji,jj,jk)                                    & 
     1334                    &                * ( 1._wp + ssh(ji,jj,Kmm_a)                       & 
     1335                    &                          * ssmask(ji,jj)                          & 
     1336                    &                          / ( ht_0(ji,jj)-1._wp + ssmask(ji,jj) ) ) 
    13231337               END DO 
    13241338            END DO 
     
    13301344         ! Save "old" scale factor (prior update) for subsequent asselin correction 
    13311345         ! of prognostic variables 
    1332          e3t_a(i1:i2,j1:j2,1:jpkm1) = e3t_n(i1:i2,j1:j2,1:jpkm1) 
    1333  
    1334          ! One should also save e3t_b, but lacking of workspace... 
    1335 !         hdivn(i1:i2,j1:j2,1:jpkm1)   = e3t_b(i1:i2,j1:j2,1:jpkm1) 
     1346         e3t(i1:i2,j1:j2,1:jpkm1,Krhs_a) = e3t(i1:i2,j1:j2,1:jpkm1,Kmm_a) 
     1347 
     1348         ! One should also save e3t(:,:,:,Kbb_a), but lacking of workspace... 
     1349!         hdiv(i1:i2,j1:j2,1:jpkm1)   = e3t(i1:i2,j1:j2,1:jpkm1,Kbb_a) 
    13361350 
    13371351         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN 
     
    13391353               DO jj=j1,j2 
    13401354                  DO ji=i1,i2 
    1341                      e3t_b(ji,jj,jk) =  e3t_b(ji,jj,jk) & 
    1342                            & + atfp * ( ptab(ji,jj,jk) - e3t_n(ji,jj,jk) ) 
     1355                     e3t(ji,jj,jk,Kbb_a) =             e3t(ji,jj,jk,Kbb_a)                  & 
     1356                       &                   + atfp * ( ptab(ji,jj,jk) - e3t(ji,jj,jk,Kmm_a) ) 
    13431357                  END DO 
    13441358               END DO 
    13451359            END DO 
    13461360            ! 
    1347             e3w_b  (i1:i2,j1:j2,1) = e3w_0(i1:i2,j1:j2,1) + e3t_b(i1:i2,j1:j2,1) - e3t_0(i1:i2,j1:j2,1) 
    1348             gdepw_b(i1:i2,j1:j2,1) = 0.0_wp 
    1349             gdept_b(i1:i2,j1:j2,1) = 0.5_wp * e3w_b(i1:i2,j1:j2,1) 
     1361            e3w  (i1:i2,j1:j2,1,Kbb_a) = e3w_0(i1:i2,j1:j2,1) + e3t(i1:i2,j1:j2,1,Kbb_a) - e3t_0(i1:i2,j1:j2,1) 
     1362            gdepw(i1:i2,j1:j2,1,Kbb_a) = 0.0_wp 
     1363            gdept(i1:i2,j1:j2,1,Kbb_a) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kbb_a) 
    13501364            ! 
    13511365            DO jk = 2, jpk 
     
    13531367                  DO ji = i1,i2             
    13541368                     zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
    1355                      e3w_b(ji,jj,jk)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) *        &  
    1356                      &                                        ( e3t_b(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) )  & 
    1357                      &                                  +            0.5_wp * tmask(ji,jj,jk)   *        & 
    1358                      &                                        ( e3t_b(ji,jj,jk  ) - e3t_0(ji,jj,jk  ) ) 
    1359                      gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1) 
    1360                      gdept_b(ji,jj,jk) =      zcoef  * ( gdepw_b(ji,jj,jk  ) + 0.5 * e3w_b(ji,jj,jk))  & 
    1361                          &               + (1-zcoef) * ( gdept_b(ji,jj,jk-1) +       e3w_b(ji,jj,jk))  
     1369                     e3w(ji,jj,jk,Kbb_a)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) )          &  
     1370                       &                                    * ( e3t(ji,jj,jk-1,Kbb_a) - e3t_0(ji,jj,jk-1) )  & 
     1371                       &                                    + 0.5_wp * tmask(ji,jj,jk)                       & 
     1372                       &                                    * ( e3t(ji,jj,jk  ,Kbb_a) - e3t_0(ji,jj,jk  ) ) 
     1373                     gdepw(ji,jj,jk,Kbb_a) =                       gdepw(ji,jj,jk-1,Kbb_a) +       e3t(ji,jj,jk-1,Kbb_a) 
     1374                     gdept(ji,jj,jk,Kbb_a) =            zcoef  * ( gdepw(ji,jj,jk  ,Kbb_a) + 0.5 * e3w(ji,jj,jk  ,Kbb_a) )  & 
     1375                       &                     + (1._wp - zcoef) * ( gdept(ji,jj,jk-1,Kbb_a) +       e3w(ji,jj,jk  ,Kbb_a) )  
    13621376                  END DO 
    13631377               END DO 
     
    13701384         ! 
    13711385         ! Update vertical scale factor at T-points: 
    1372          e3t_n(i1:i2,j1:j2,1:jpkm1) = ptab(i1:i2,j1:j2,1:jpkm1) 
     1386         e3t(i1:i2,j1:j2,1:jpkm1,Kmm_a) = ptab(i1:i2,j1:j2,1:jpkm1) 
    13731387         ! 
    13741388         ! Update total depth: 
    1375          ht_n(i1:i2,j1:j2) = 0._wp 
     1389         ht(i1:i2,j1:j2) = 0._wp 
    13761390         DO jk = 1, jpkm1 
    1377             ht_n(i1:i2,j1:j2) = ht_n(i1:i2,j1:j2) + e3t_n(i1:i2,j1:j2,jk) * tmask(i1:i2,j1:j2,jk) 
     1391            ht(i1:i2,j1:j2) = ht(i1:i2,j1:j2) + e3t(i1:i2,j1:j2,jk,Kmm_a) * tmask(i1:i2,j1:j2,jk) 
    13781392         END DO 
    13791393         ! 
    13801394         ! Update vertical scale factor at W-points and depths: 
    1381          e3w_n (i1:i2,j1:j2,1) = e3w_0(i1:i2,j1:j2,1) + e3t_n(i1:i2,j1:j2,1) - e3t_0(i1:i2,j1:j2,1) 
    1382          gdept_n(i1:i2,j1:j2,1) = 0.5_wp * e3w_n(i1:i2,j1:j2,1) 
    1383          gdepw_n(i1:i2,j1:j2,1) = 0.0_wp 
    1384          gde3w_n(i1:i2,j1:j2,1) = gdept_n(i1:i2,j1:j2,1) - (ht_n(i1:i2,j1:j2)-ht_0(i1:i2,j1:j2)) ! Last term in the rhs is ssh 
     1395           e3w(i1:i2,j1:j2,1,Kmm_a) = e3w_0(i1:i2,j1:j2,1) + e3t(i1:i2,j1:j2,1,Kmm_a) - e3t_0(i1:i2,j1:j2,1) 
     1396         gdept(i1:i2,j1:j2,1,Kmm_a) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kmm_a) 
     1397         gdepw(i1:i2,j1:j2,1,Kmm_a) = 0.0_wp 
     1398         gde3w(i1:i2,j1:j2,1      ) = gdept(i1:i2,j1:j2,1,Kmm_a) - ( ht(i1:i2,j1:j2) - ht_0(i1:i2,j1:j2) ) ! Last term in the rhs is ssh 
    13851399         ! 
    13861400         DO jk = 2, jpk 
     
    13881402               DO ji = i1,i2             
    13891403               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
    1390                e3w_n(ji,jj,jk)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * ( e3t_n(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) )   & 
    1391                &                                  +            0.5_wp * tmask(ji,jj,jk)   * ( e3t_n(ji,jj,jk  ) - e3t_0(ji,jj,jk  ) ) 
    1392                gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 
    1393                gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  & 
    1394                    &               + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk))  
    1395                gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - (ht_n(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh 
     1404                e3w(ji,jj,jk,Kmm_a)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) )          & 
     1405                 &                                     * ( e3t(ji,jj,jk-1,Kmm_a) - e3t_0(ji,jj,jk-1) )  & 
     1406                 &                                     + 0.5_wp * tmask(ji,jj,jk)                       & 
     1407                 &                                     * ( e3t(ji,jj,jk  ,Kmm_a) - e3t_0(ji,jj,jk  ) ) 
     1408               gdepw(ji,jj,jk,Kmm_a) =                         gdepw(ji,jj,jk-1,Kmm_a) +       e3t(ji,jj,jk-1,Kmm_a) 
     1409               gdept(ji,jj,jk,Kmm_a) =             zcoef   * ( gdepw(ji,jj,jk  ,Kmm_a) + 0.5 * e3w(ji,jj,jk  ,Kmm_a) )  & 
     1410                 &                     + ( 1._wp - zcoef ) * ( gdept(ji,jj,jk-1,Kmm_a) +       e3w(ji,jj,jk  ,Kmm_a) )  
     1411               gde3w(ji,jj,jk      ) =                         gdept(ji,jj,jk  ,Kmm_a) - ( ht(ji,jj)-ht_0(ji,jj) )    ! Last term in the rhs is ssh 
    13961412               END DO 
    13971413            END DO 
     
    13991415         ! 
    14001416         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
    1401             e3t_b (i1:i2,j1:j2,1:jpk)  = e3t_n (i1:i2,j1:j2,1:jpk) 
    1402             e3w_b (i1:i2,j1:j2,1:jpk)  = e3w_n (i1:i2,j1:j2,1:jpk) 
    1403             gdepw_b(i1:i2,j1:j2,1:jpk) = gdepw_n(i1:i2,j1:j2,1:jpk) 
    1404             gdept_b(i1:i2,j1:j2,1:jpk) = gdept_n(i1:i2,j1:j2,1:jpk) 
     1417              e3t(i1:i2,j1:j2,1:jpk,Kbb_a) =   e3t(i1:i2,j1:j2,1:jpk,Kmm_a) 
     1418              e3w(i1:i2,j1:j2,1:jpk,Kbb_a) =   e3w(i1:i2,j1:j2,1:jpk,Kmm_a) 
     1419            gdepw(i1:i2,j1:j2,1:jpk,Kbb_a) = gdepw(i1:i2,j1:j2,1:jpk,Kmm_a) 
     1420            gdept(i1:i2,j1:j2,1:jpk,Kbb_a) = gdept(i1:i2,j1:j2,1:jpk,Kmm_a) 
    14051421         ENDIF 
    14061422         ! 
Note: See TracChangeset for help on using the changeset viewer.