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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 9031 for branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90 – NEMO

Ignore:
Timestamp:
2017-12-14T11:10:02+01:00 (6 years ago)
Author:
timgraham
Message:

Resolved AGRIF conflicts

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90

    r9019 r9031  
    11#define TWO_WAY        /* TWO WAY NESTING */ 
    2 #undef DECAL_FEEDBACK /* SEPARATION of INTERFACES*/ 
     2#define DECAL_FEEDBACK /* SEPARATION of INTERFACES*/ 
     3#undef VOL_REFLUX      /* VOLUME REFLUXING*/ 
    34  
    45MODULE agrif_opa_update 
    … …  
    2324   USE in_out_manager ! I/O manager 
    2425   USE lib_mpp        ! MPP library 
     26   USE domvvl         ! Need interpolation routines  
    2527 
    2628   IMPLICIT NONE 
    … …  
    3638CONTAINS 
    3739 
    38    RECURSIVE SUBROUTINE Agrif_Update_Tra( ) 
     40   SUBROUTINE Agrif_Update_Tra( ) 
    3941      !!---------------------------------------------------------------------- 
    4042      !!                   *** ROUTINE Agrif_Update_Tra *** 
    … …  
    4446      ! 
    4547#if defined TWO_WAY   
    46       IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update tracers  from grid Number',Agrif_Fixed(), 'nbcline', nbcline 
     48      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update tracers  from grid Number',Agrif_Fixed() 
    4749 
    4850      Agrif_UseSpecialValueInUpdate = .TRUE. 
    4951      Agrif_SpecialValueFineGrid    = 0._wp 
    5052      !  
    51       IF (MOD(nbcline,nbclineupdate) == 0) THEN 
    5253# if ! defined DECAL_FEEDBACK 
    53          CALL Agrif_Update_Variable(tsn_id, procname=updateTS) 
     54      CALL Agrif_Update_Variable(tsn_id, procname=updateTS) 
     55! near boundary update: 
     56!      CALL Agrif_Update_Variable(tsn_id,locupdate=(/0,2/), procname=updateTS) 
    5457# else 
    55          CALL Agrif_Update_Variable(tsn_id, locupdate=(/1,0/),procname=updateTS) 
     58      CALL Agrif_Update_Variable(tsn_id, locupdate=(/1,0/),procname=updateTS) 
     59! near boundary update: 
     60!      CALL Agrif_Update_Variable(tsn_id,locupdate=(/1,2/), procname=updateTS) 
    5661# endif 
    57       ELSE 
    58 # if ! defined DECAL_FEEDBACK 
    59          CALL Agrif_Update_Variable(tsn_id,locupdate=(/0,2/), procname=updateTS) 
    60 # else 
    61          CALL Agrif_Update_Variable(tsn_id,locupdate=(/1,2/), procname=updateTS) 
    62 # endif 
    63       ENDIF 
    6462      ! 
    6563      Agrif_UseSpecialValueInUpdate = .FALSE. 
    6664      ! 
    67       IF ( lk_agrif_doupd ) THEN ! Initialisation: do recursive update: 
    68          CALL Agrif_ChildGrid_To_ParentGrid() 
    69          CALL Agrif_Update_Tra() 
    70          CALL Agrif_ParentGrid_To_ChildGrid() 
    71       ENDIF 
    72       ! 
    7365#endif 
    7466      ! 
    7567   END SUBROUTINE Agrif_Update_Tra 
    7668 
    77  
    78    RECURSIVE SUBROUTINE Agrif_Update_Dyn( ) 
     69   SUBROUTINE Agrif_Update_Dyn( ) 
    7970      !!---------------------------------------------------------------------- 
    8071      !!                   *** ROUTINE Agrif_Update_Dyn *** 
    … …  
    8475      ! 
    8576#if defined TWO_WAY 
    86       IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update momentum from grid Number',Agrif_Fixed(), 'nbcline', nbcline 
     77      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update momentum from grid Number',Agrif_Fixed() 
    8778 
    8879      Agrif_UseSpecialValueInUpdate = .FALSE. 
    8980      Agrif_SpecialValueFineGrid = 0. 
    9081      !      
    91       IF (mod(nbcline,nbclineupdate) == 0) THEN 
    9282# if ! defined DECAL_FEEDBACK 
    93          CALL Agrif_Update_Variable(un_update_id,procname = updateU) 
    94          CALL Agrif_Update_Variable(vn_update_id,procname = updateV) 
     83      CALL Agrif_Update_Variable(un_update_id,procname = updateU) 
     84      CALL Agrif_Update_Variable(vn_update_id,procname = updateV) 
     85! near boundary update: 
     86!      CALL Agrif_Update_Variable(un_update_id,locupdate=(/0,1/),procname = updateU) 
     87!      CALL Agrif_Update_Variable(vn_update_id,locupdate=(/0,1/),procname = updateV) 
    9588# else 
    96          CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU) 
    97          CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV) 
     89      CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU) 
     90      CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV) 
     91! near boundary update: 
     92!      CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,1/),locupdate2=(/1,1/),procname = updateU) 
     93!      CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updateV) 
    9894# endif 
    99       ELSE 
    100 # if ! defined DECAL_FEEDBACK 
    101          CALL Agrif_Update_Variable(un_update_id,locupdate=(/0,1/),procname = updateU) 
    102          CALL Agrif_Update_Variable(vn_update_id,locupdate=(/0,1/),procname = updateV)          
    103 # else 
    104          CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,1/),locupdate2=(/1,1/),procname = updateU) 
    105          CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updateV) 
    106 # endif 
    107       ENDIF 
    10895 
    10996# if ! defined DECAL_FEEDBACK 
    … …  
    114101      CALL Agrif_Update_Variable(e2v_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV2d)   
    115102# endif 
    116  
     103      ! 
     104# if ! defined DECAL_FEEDBACK 
     105      ! Account for updated thicknesses at boundary edges 
     106      IF (.NOT.ln_linssh) THEN 
     107!         CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_u_bdy) 
     108!         CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_v_bdy) 
     109      ENDIF 
     110# endif 
     111      !  
    117112      IF ( ln_dynspg_ts .AND. ln_bt_fw ) THEN 
    118113         ! Update time integrated transports 
    119          IF (mod(nbcline,nbclineupdate) == 0) THEN 
    120114#  if ! defined DECAL_FEEDBACK 
    121             CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b) 
    122             CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b) 
     115         CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b) 
     116         CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b) 
    123117#  else 
    124             CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateub2b) 
    125             CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updatevb2b) 
     118         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateub2b) 
     119         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updatevb2b) 
    126120#  endif 
    127          ELSE 
    128 #  if ! defined DECAL_FEEDBACK 
    129             CALL Agrif_Update_Variable(ub2b_update_id,locupdate=(/0,1/),procname = updateub2b) 
    130             CALL Agrif_Update_Variable(vb2b_update_id,locupdate=(/0,1/),procname = updatevb2b) 
    131 #  else 
    132             CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,1/),locupdate2=(/1,1/),procname = updateub2b) 
    133             CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updatevb2b) 
    134 #  endif 
    135          ENDIF 
    136121      END IF 
    137       ! 
    138       nbcline = nbcline + 1 
     122#endif 
     123      ! 
     124   END SUBROUTINE Agrif_Update_Dyn 
     125 
     126   SUBROUTINE Agrif_Update_ssh( ) 
     127      !!--------------------------------------------- 
     128      !!   *** ROUTINE Agrif_Update_ssh *** 
     129      !!--------------------------------------------- 
     130      !  
     131      IF (Agrif_Root()) RETURN 
     132      ! 
     133#if defined TWO_WAY 
    139134      ! 
    140135      Agrif_UseSpecialValueInUpdate = .TRUE. 
    … …  
    145140      CALL Agrif_Update_Variable(sshn_id,locupdate=(/1,0/),procname = updateSSH) 
    146141# endif 
     142      ! 
    147143      Agrif_UseSpecialValueInUpdate = .FALSE. 
    148       !  
     144      ! 
     145#  if defined DECAL_FEEDBACK && defined VOL_REFLUX 
     146      IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN 
     147         ! Refluxing on ssh: 
     148         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0, 0/),locupdate2=(/1, 1/),procname = reflux_sshu) 
     149         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1, 1/),locupdate2=(/0, 0/),procname = reflux_sshv) 
     150      END IF 
     151#  endif 
     152      ! 
    149153#endif 
    150154      ! 
    151       ! Do recursive update: 
    152       IF ( lk_agrif_doupd ) THEN ! Initialisation: do recursive update: 
    153          CALL Agrif_ChildGrid_To_ParentGrid() 
    154          CALL Agrif_Update_Dyn() 
    155          CALL Agrif_ParentGrid_To_ChildGrid() 
    156       ENDIF 
    157       ! 
    158    END SUBROUTINE Agrif_Update_Dyn 
    159  
     155   END SUBROUTINE Agrif_Update_ssh 
     156 
     157 
     158   SUBROUTINE Agrif_Update_Tke( kt ) 
     159      !!--------------------------------------------- 
     160      !!   *** ROUTINE Agrif_Update_Tke *** 
     161      !!--------------------------------------------- 
     162      !! 
     163      INTEGER, INTENT(in) :: kt  
     164      !  
     165      IF (Agrif_Root()) RETURN 
     166      !        
     167      IF( ( Agrif_NbStepint() .NE. 0 ) .AND. (kt /= 0) ) RETURN 
     168#  if defined TWO_WAY 
     169 
     170      Agrif_UseSpecialValueInUpdate = .TRUE. 
     171      Agrif_SpecialValueFineGrid = 0. 
     172 
     173      CALL Agrif_Update_Variable( en_id, locupdate=(/0,0/), procname=updateEN  ) 
     174      CALL Agrif_Update_Variable(avt_id, locupdate=(/0,0/), procname=updateAVT ) 
     175      CALL Agrif_Update_Variable(avm_id, locupdate=(/0,0/), procname=updateAVM ) 
     176 
     177      Agrif_UseSpecialValueInUpdate = .FALSE. 
     178 
     179#  endif 
     180       
     181   END SUBROUTINE Agrif_Update_Tke 
     182 
     183 
     184   SUBROUTINE Agrif_Update_vvl( ) 
     185      !!--------------------------------------------- 
     186      !!   *** ROUTINE Agrif_Update_vvl *** 
     187      !!--------------------------------------------- 
     188      ! 
     189      IF (Agrif_Root()) RETURN 
     190      ! 
     191#if defined TWO_WAY   
     192      ! 
     193      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update e3 from grid Number',Agrif_Fixed(), 'Step', Agrif_Nb_Step() 
     194      ! 
     195      Agrif_UseSpecialValueInUpdate = .TRUE. 
     196      Agrif_SpecialValueFineGrid = 0. 
     197      !  
     198      ! No interface separation here, update vertical grid at T points  
     199      ! everywhere over the overlapping regions (one account for refluxing in that case): 
     200      CALL Agrif_Update_Variable(e3t_id, procname=updatee3t)  
     201      ! 
     202      Agrif_UseSpecialValueInUpdate = .FALSE. 
     203      ! 
     204      CALL Agrif_ChildGrid_To_ParentGrid() 
     205      CALL dom_vvl_update_UVF 
     206      CALL Agrif_ParentGrid_To_ChildGrid() 
     207      ! 
     208#endif 
     209      ! 
     210   END SUBROUTINE Agrif_Update_vvl 
     211 
     212   SUBROUTINE dom_vvl_update_UVF 
     213      !!--------------------------------------------- 
     214      !!       *** ROUTINE dom_vvl_update_UVF *** 
     215      !!--------------------------------------------- 
     216      !! 
     217      INTEGER :: jk 
     218      REAL(wp):: zcoef 
     219      !!--------------------------------------------- 
     220 
     221      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Finalize e3 on grid Number', & 
     222                  & Agrif_Fixed(), 'Step', Agrif_Nb_Step() 
     223 
     224      ! Save "old" scale factor (prior update) for subsequent asselin correction 
     225      ! of prognostic variables 
     226      ! ----------------------- 
     227      ! 
     228      e3u_a(:,:,:) = e3u_n(:,:,:) 
     229      e3v_a(:,:,:) = e3v_n(:,:,:) 
     230!      ua(:,:,:) = e3u_b(:,:,:) 
     231!      va(:,:,:) = e3v_b(:,:,:) 
     232      hu_a(:,:) = hu_n(:,:) 
     233      hv_a(:,:) = hv_n(:,:) 
     234 
     235      ! 1) NOW fields 
     236      !-------------- 
     237       
     238         ! Vertical scale factor interpolations 
     239         ! ------------------------------------ 
     240      CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:) ,  'U' ) 
     241      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:) ,  'V' ) 
     242      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:) ,  'F' ) 
     243 
     244      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 
     245      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 
     246 
     247         ! Update total depths: 
     248         ! -------------------- 
     249      hu_n(:,:) = 0._wp                        ! Ocean depth at U-points 
     250      hv_n(:,:) = 0._wp                        ! Ocean depth at V-points 
     251      DO jk = 1, jpkm1 
     252         hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk) 
     253         hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk) 
     254      END DO 
     255      !                                        ! Inverse of the local depth 
     256      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) 
     257      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) 
     258 
     259 
     260      ! 2) BEFORE fields: 
     261      !------------------ 
     262!      IF (     (.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_exp)) & 
     263!         & .OR.(.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_ts    & 
     264!         & .AND.(.NOT.ln_bt_fw)))) THEN 
     265      IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN 
     266         ! 
     267         ! Vertical scale factor interpolations 
     268         ! ------------------------------------ 
     269         CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:),  'U'  ) 
     270         CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:),  'V'  ) 
     271 
     272         CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 
     273         CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 
     274 
     275         ! Update total depths: 
     276         ! -------------------- 
     277         hu_b(:,:) = 0._wp                     ! Ocean depth at U-points 
     278         hv_b(:,:) = 0._wp                     ! Ocean depth at V-points 
     279         DO jk = 1, jpkm1 
     280            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk) 
     281            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk) 
     282         END DO 
     283         !                                     ! Inverse of the local depth 
     284         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) 
     285         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) 
     286      ENDIF 
     287      ! 
     288   END SUBROUTINE dom_vvl_update_UVF 
     289 
     290#if defined key_vertical 
    160291 
    161292   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    162293      !!---------------------------------------------------------------------- 
    163294      !!           *** ROUTINE updateT *** 
    164       !!---------------------------------------------------------------------- 
    165       INTEGER                                    , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2 
    166       REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   tabres 
    167       LOGICAL                                    , INTENT(in   ) ::   before 
    168       ! 
    169       INTEGER :: ji, jj, jk, jn 
    170       !!---------------------------------------------------------------------- 
    171       ! 
    172       IF( before ) THEN 
    173          DO jn = n1, n2 
    174             DO jk = k1, k2 
    175                DO jj = j1, j2 
    176                   DO ji = i1, i2 
    177                      tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) 
     295      !!--------------------------------------------- 
     296      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2 
     297      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
     298      LOGICAL, INTENT(in) :: before 
     299      !! 
     300      INTEGER :: ji,jj,jk,jn 
     301      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: tabres_child 
     302      REAL(wp) :: h_in(k1:k2) 
     303      REAL(wp) :: h_out(1:jpk) 
     304      INTEGER  :: N_in, N_out 
     305      REAL(wp) :: h_diff 
     306      REAL(wp) :: zrho_xy 
     307      REAL(wp) :: tabin(k1:k2,n1:n2) 
     308      !!--------------------------------------------- 
     309      ! 
     310      IF (before) THEN 
     311         AGRIF_SpecialValue = -999._wp 
     312         zrho_xy = Agrif_rhox() * Agrif_rhoy()  
     313         DO jn = n1,n2-1 
     314            DO jk=k1,k2 
     315               DO jj=j1,j2 
     316                  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 
    178319                  END DO 
    179320               END DO 
    180321            END DO 
    181322         END DO 
     323         DO jk=k1,k2 
     324            DO jj=j1,j2 
     325               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 
     328               END DO 
     329            END DO 
     330         END DO 
    182331      ELSE 
     332         tabres_child(:,:,:,:) = 0. 
     333         AGRIF_SpecialValue = 0._wp 
     334         DO jj=j1,j2 
     335            DO ji=i1,i2 
     336               N_in = 0 
     337               DO jk=k1,k2 !k2 = jpk of child grid 
     338                  IF (tabres(ji,jj,jk,n2) == 0  ) EXIT 
     339                  N_in = N_in + 1 
     340                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 
     341                  h_in(N_in) = tabres(ji,jj,jk,n2) 
     342               ENDDO 
     343               N_out = 0 
     344               DO jk=1,jpk ! jpk of parent grid 
     345                  IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF 
     346                  N_out = N_out + 1 
     347                  h_out(N_out) = e3t_n(ji,jj,jk)  
     348               ENDDO 
     349               IF (N_in > 0) THEN !Remove this? 
     350                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     351                  IF (h_diff < -1.e-4) THEN 
     352                     print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,sum(h_in(1:N_in)),sum(h_out(1:N_out)) 
     353                     print *,h_in(1:N_in) 
     354                     print *,h_out(1:N_out) 
     355                     STOP 
     356                  ENDIF 
     357                  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) 
     359                  ENDDO 
     360               ENDIF 
     361            ENDDO 
     362         ENDDO 
     363 
    183364         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
    184365            ! Add asselin part 
    185             DO jn = n1,n2 
     366            DO jn = n1,n2-1 
     367               DO jk=1,jpk 
     368                  DO jj=j1,j2 
     369                     DO ji=i1,i2 
     370                        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) 
     374                        ENDIF 
     375                     ENDDO 
     376                  ENDDO 
     377               ENDDO 
     378            ENDDO 
     379         ENDIF 
     380         DO jn = n1,n2-1 
     381            DO jk=1,jpk 
     382               DO jj=j1,j2 
     383                  DO ji=i1,i2 
     384                     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) 
     386                     END IF 
     387                  END DO 
     388               END DO 
     389            END DO 
     390         END DO 
     391      ENDIF 
     392      !  
     393   END SUBROUTINE updateTS 
     394 
     395# else 
     396 
     397   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
     398      !!--------------------------------------------- 
     399      !!           *** ROUTINE updateT *** 
     400      !!--------------------------------------------- 
     401      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2 
     402      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
     403      LOGICAL, INTENT(in) :: before 
     404      !! 
     405      INTEGER :: ji,jj,jk,jn 
     406      REAL(wp) :: ztb, ztnu, ztno 
     407      !!--------------------------------------------- 
     408      ! 
     409      IF (before) THEN 
     410         DO jn = 1,jpts 
     411            DO jk=k1,k2 
     412               DO jj=j1,j2 
     413                  DO ji=i1,i2 
     414!> 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!< jc tmp 
     418                  END DO 
     419               END DO 
     420            END DO 
     421         END DO 
     422      ELSE 
     423!> jc tmp 
     424         DO jn = 1,jpts 
     425            tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) & 
     426                                         & * tmask(i1:i2,j1:j2,k1:k2) 
     427         ENDDO 
     428!< jc tmp 
     429         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
     430            ! Add asselin part 
     431            DO jn = 1,jpts 
    186432               DO jk = k1, k2 
    187433                  DO jj = j1, j2 
    188434                     DO ji = i1, i2 
    189435                        IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN 
    190                            tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) &  
    191                               &             + atfp * ( tabres(ji,jj,jk,jn) - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 
     436                           ztb  = tsb(ji,jj,jk,jn) * e3t_b(ji,jj,jk) ! fse3t_b prior update should be used 
     437                           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) 
    192441                        ENDIF 
    193442                     END DO 
    … …  
    196445            END DO 
    197446         ENDIF 
    198          DO jn = n1,n2 
     447         DO jn = 1,jpts 
    199448            DO jk=k1,k2 
    200449               DO jj=j1,j2 
    201450                  DO ji=i1,i2 
    202451                     IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN  
    203                         tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) * tmask(ji,jj,jk) 
     452                        tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) / e3t_n(ji,jj,jk) 
    204453                     END IF 
    205454                  END DO 
    … …  
    207456            END DO 
    208457         END DO 
     458         ! 
     459         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) 
     461         ENDIF 
     462         ! 
    209463      ENDIF 
    210464      !  
    211465   END SUBROUTINE updateTS 
    212466 
    213  
    214    SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, before ) 
     467# endif 
     468 
     469# if defined key_vertical 
     470 
     471   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    215472      !!--------------------------------------------- 
    216473      !!           *** ROUTINE updateu *** 
    217474      !!--------------------------------------------- 
    218       INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2 
    219       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres 
    220       LOGICAL                               , INTENT(in   ) :: before 
     475      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
     476      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
     477      LOGICAL                                     , INTENT(in   ) :: before 
    221478      ! 
    222479      INTEGER ::   ji, jj, jk 
    223480      REAL(wp)::   zrhoy 
     481! VERTICAL REFINEMENT BEGIN 
     482      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 
     483      REAL(wp) :: h_in(k1:k2) 
     484      REAL(wp) :: h_out(1:jpk) 
     485      INTEGER  :: N_in, N_out 
     486      REAL(wp) :: h_diff, excess, thick 
     487      REAL(wp) :: tabin(k1:k2) 
     488! VERTICAL REFINEMENT END 
     489      !!--------------------------------------------- 
     490      !  
     491      IF( before ) THEN 
     492         zrhoy = Agrif_Rhoy() 
     493         AGRIF_SpecialValue = -999._wp 
     494         DO jk=k1,k2 
     495            DO jj=j1,j2 
     496               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 
     501               END DO 
     502            END DO 
     503         END DO 
     504      ELSE 
     505         tabres_child(:,:,:) = 0. 
     506         AGRIF_SpecialValue = 0._wp 
     507         DO jj=j1,j2 
     508            DO ji=i1,i2 
     509               N_in = 0 
     510               h_in(:) = 0._wp 
     511               tabin(:) = 0._wp 
     512               DO jk=k1,k2 !k2=jpk of child grid 
     513                  IF( tabres(ji,jj,jk,2) < -900) EXIT 
     514                  N_in = N_in + 1 
     515                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
     516                  h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj) 
     517               ENDDO 
     518               N_out = 0 
     519               DO jk=1,jpk 
     520                  IF (umask(ji,jj,jk) == 0) EXIT 
     521                  N_out = N_out + 1 
     522                  h_out(N_out) = e3u_n(ji,jj,jk) 
     523               ENDDO 
     524               IF (N_in * N_out > 0) THEN 
     525                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     526                  IF (h_diff < -1.e-4) THEN 
     527!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.  
     528!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. 
     529                     excess = 0._wp 
     530                     DO jk=N_in,1,-1 
     531                        thick = MIN(-1*h_diff, h_in(jk)) 
     532                        excess = excess + tabin(jk)*thick*e2u(ji,jj) 
     533                        tabin(jk) = tabin(jk)*(1. - thick/h_in(jk)) 
     534                        h_diff = h_diff + thick 
     535                        IF ( h_diff == 0) THEN 
     536                           N_in = jk 
     537                           h_in(jk) = h_in(jk) - thick 
     538                           EXIT 
     539                        ENDIF 
     540                     ENDDO 
     541                  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               ENDIF 
     545            ENDDO 
     546         ENDDO 
     547 
     548         DO jk=1,jpk 
     549            DO jj=j1,j2 
     550               DO ji=i1,i2 
     551                  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) 
     554                  ENDIF 
     555                  ! 
     556                  un(ji,jj,jk) = tabres_child(ji,jj,jk) * umask(ji,jj,jk) 
     557               END DO 
     558            END DO 
     559         END DO 
     560      ENDIF 
     561      !  
     562   END SUBROUTINE updateu 
     563 
     564#else 
     565 
     566   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
     567      !!--------------------------------------------- 
     568      !!           *** ROUTINE updateu *** 
     569      !!--------------------------------------------- 
     570      INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
     571      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
     572      LOGICAL                                     , INTENT(in   ) :: before 
     573      ! 
     574      INTEGER  :: ji, jj, jk 
     575      REAL(wp) :: zrhoy, zub, zunu, zuno 
    224576      !!--------------------------------------------- 
    225577      !  
    … …  
    227579         zrhoy = Agrif_Rhoy() 
    228580         DO jk = k1, k2 
    229             tabres(i1:i2,j1:j2,jk) = zrhoy * e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk) 
     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) 
    230582         END DO 
    231583      ELSE 
    … …  
    233585            DO jj=j1,j2 
    234586               DO ji=i1,i2 
    235                   tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk) 
     587                  tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e2u(ji,jj)  
    236588                  ! 
    237589                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    238                      ub(ji,jj,jk) = ub(ji,jj,jk) &  
    239                            & + atfp * ( tabres(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk) 
     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) 
     592                     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) 
    240595                  ENDIF 
    241596                  ! 
    242                   un(ji,jj,jk) = tabres(ji,jj,jk) * umask(ji,jj,jk) 
    243                END DO 
    244             END DO 
    245          END DO 
     597                  un(ji,jj,jk) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u_n(ji,jj,jk) 
     598               END DO 
     599            END DO 
     600         END DO 
     601         ! 
     602         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     603            ub(i1:i2,j1:j2,k1:k2)  = un(i1:i2,j1:j2,k1:k2) 
     604         ENDIF 
     605         ! 
    246606      ENDIF 
    247607      !  
    248608   END SUBROUTINE updateu 
    249609 
    250  
    251    SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, before ) 
    252       !!---------------------------------------------------------------------- 
    253       !!                      *** ROUTINE updatev *** 
    254       !!---------------------------------------------------------------------- 
    255       INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2 
    256       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres 
    257       LOGICAL                               , INTENT(in   ) :: before 
     610# endif 
     611 
     612   SUBROUTINE correct_u_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) 
     613      !!--------------------------------------------- 
     614      !!           *** ROUTINE correct_u_bdy *** 
     615      !!--------------------------------------------- 
     616      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
     617      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
     618      LOGICAL                                     , INTENT(in   ) :: before 
     619      INTEGER                                     , INTENT(in)    :: nb, ndir 
    258620      !! 
     621      LOGICAL :: western_side, eastern_side  
     622      ! 
     623      INTEGER  :: jj, jk 
     624      REAL(wp) :: zcor 
     625      !!--------------------------------------------- 
     626      !  
     627      IF( .NOT.before ) THEN 
     628         ! 
     629         western_side  = (nb == 1).AND.(ndir == 1) 
     630         eastern_side  = (nb == 1).AND.(ndir == 2) 
     631         ! 
     632         IF (western_side) THEN 
     633            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 
     636               DO jk=1,jpkm1 
     637                  un(i1-1,jj,jk) = un(i1-1,jj,jk) + zcor * umask(i1-1,jj,jk) 
     638               END DO  
     639            END DO 
     640         ENDIF 
     641         ! 
     642         IF (eastern_side) THEN 
     643            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 
     646               DO jk=1,jpkm1 
     647                  un(i2+1,jj,jk) = un(i2+1,jj,jk) + zcor * umask(i2+1,jj,jk) 
     648               END DO  
     649            END DO 
     650         ENDIF 
     651         ! 
     652      ENDIF 
     653      !  
     654   END SUBROUTINE correct_u_bdy 
     655 
     656# if  defined key_vertical 
     657 
     658   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
     659      !!--------------------------------------------- 
     660      !!           *** ROUTINE updatev *** 
     661      !!--------------------------------------------- 
     662      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
     663      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
     664      LOGICAL                                     , INTENT(in   ) :: before 
     665      ! 
    259666      INTEGER  ::   ji, jj, jk 
    260667      REAL(wp) ::   zrhox 
    261       !!---------------------------------------------------------------------- 
     668! VERTICAL REFINEMENT BEGIN 
     669      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 
     670      REAL(wp) :: h_in(k1:k2) 
     671      REAL(wp) :: h_out(1:jpk) 
     672      INTEGER :: N_in, N_out 
     673      REAL(wp) :: h_diff, excess, thick 
     674      REAL(wp) :: tabin(k1:k2) 
     675! VERTICAL REFINEMENT END 
     676      !!---------------------------------------------       
    262677      ! 
    263678      IF( before ) THEN 
    264679         zrhox = Agrif_Rhox() 
    265          DO jk = k1, k2 
    266             DO jj = j1, j2 
    267                DO ji = i1, i2 
    268                   tabres(ji,jj,jk) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk) 
     680         AGRIF_SpecialValue = -999._wp 
     681         DO jk=k1,k2 
     682            DO jj=j1,j2 
     683               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) & 
     687                                       + (vmask(ji,jj,jk)-1)*999._wp 
    269688               END DO 
    270689            END DO 
    271690         END DO 
    272691      ELSE 
    273          DO jk = k1, k2 
    274             DO jj = j1, j2 
    275                DO ji = i1, i2 
    276                   tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk) 
     692         tabres_child(:,:,:) = 0. 
     693         AGRIF_SpecialValue = 0._wp 
     694         DO jj=j1,j2 
     695            DO ji=i1,i2 
     696               N_in = 0 
     697               DO jk=k1,k2 
     698                  IF (tabres(ji,jj,jk,2) < -900) EXIT 
     699                  N_in = N_in + 1 
     700                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
     701                  h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj) 
     702               ENDDO 
     703               N_out = 0 
     704               DO jk=1,jpk 
     705                  IF (vmask(ji,jj,jk) == 0) EXIT 
     706                  N_out = N_out + 1 
     707                  h_out(N_out) = e3v_n(ji,jj,jk) 
     708               ENDDO 
     709               IF (N_in * N_out > 0) THEN 
     710                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     711                  IF (h_diff < -1.e-4) then 
     712!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.  
     713!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. 
     714                     excess = 0._wp 
     715                     DO jk=N_in,1,-1 
     716                        thick = MIN(-1*h_diff, h_in(jk)) 
     717                        excess = excess + tabin(jk)*thick*e2u(ji,jj) 
     718                        tabin(jk) = tabin(jk)*(1. - thick/h_in(jk)) 
     719                        h_diff = h_diff + thick 
     720                        IF ( h_diff == 0) THEN 
     721                           N_in = jk 
     722                           h_in(jk) = h_in(jk) - thick 
     723                           EXIT 
     724                        ENDIF 
     725                     ENDDO 
     726                  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) 
     728                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out)) 
     729               ENDIF 
     730            ENDDO 
     731         ENDDO 
     732 
     733         DO jk=1,jpk 
     734            DO jj=j1,j2 
     735               DO ji=i1,i2 
    277736                  ! 
    278737                  IF( .NOT.(lk_agrif_fstep.AND.(neuler==0)) ) THEN ! Add asselin part 
    279738                     vb(ji,jj,jk) = vb(ji,jj,jk) &  
    280                         &         + atfp * ( tabres(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk) 
     739                           & + atfp * ( tabres_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk) 
    281740                  ENDIF 
    282741                  ! 
    283                   vn(ji,jj,jk) = tabres(ji,jj,jk) * vmask(ji,jj,jk) 
     742                  vn(ji,jj,jk) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk) 
    284743               END DO 
    285744            END DO 
    … …  
    288747      !  
    289748   END SUBROUTINE updatev 
     749 
     750# else 
     751 
     752   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before) 
     753      !!--------------------------------------------- 
     754      !!           *** ROUTINE updatev *** 
     755      !!--------------------------------------------- 
     756      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
     757      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
     758      LOGICAL                                     , INTENT(in   ) :: before 
     759      ! 
     760      INTEGER  :: ji, jj, jk 
     761      REAL(wp) :: zrhox, zvb, zvnu, zvno 
     762      !!---------------------------------------------       
     763      ! 
     764      IF (before) THEN 
     765         zrhox = Agrif_Rhox() 
     766         DO jk=k1,k2 
     767            DO jj=j1,j2 
     768               DO ji=i1,i2 
     769                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk) 
     770               END DO 
     771            END DO 
     772         END DO 
     773      ELSE 
     774         DO jk=k1,k2 
     775            DO jj=j1,j2 
     776               DO ji=i1,i2 
     777                  tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e1v(ji,jj) 
     778                  ! 
     779                  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) 
     782                     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) 
     785                  ENDIF 
     786                  ! 
     787                  vn(ji,jj,jk) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v_n(ji,jj,jk) 
     788               END DO 
     789            END DO 
     790         END DO 
     791         ! 
     792         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     793            vb(i1:i2,j1:j2,k1:k2)  = vn(i1:i2,j1:j2,k1:k2) 
     794         ENDIF 
     795         ! 
     796      ENDIF 
     797      !  
     798   END SUBROUTINE updatev 
     799 
     800# endif 
     801 
     802   SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) 
     803      !!--------------------------------------------- 
     804      !!           *** ROUTINE correct_u_bdy *** 
     805      !!--------------------------------------------- 
     806      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
     807      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
     808      LOGICAL                                     , INTENT(in   ) :: before 
     809      INTEGER                                     , INTENT(in)    :: nb, ndir 
     810      !! 
     811      LOGICAL :: southern_side, northern_side  
     812      ! 
     813      INTEGER  :: ji, jk 
     814      REAL(wp) :: zcor 
     815      !!--------------------------------------------- 
     816      !  
     817      IF( .NOT.before ) THEN 
     818         ! 
     819         southern_side = (nb == 2).AND.(ndir == 1) 
     820         northern_side = (nb == 2).AND.(ndir == 2) 
     821         ! 
     822         IF (southern_side) THEN 
     823            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 
     826               DO jk=1,jpkm1 
     827                  vn(ji,j1-1,jk) = vn(ji,j1-1,jk) + zcor * vmask(ji,j1-1,jk) 
     828               END DO  
     829            END DO 
     830         ENDIF 
     831         ! 
     832         IF (northern_side) THEN 
     833            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 
     836               DO jk=1,jpkm1 
     837                  vn(ji,j2+1,jk) = vn(ji,j2+1,jk) + zcor * vmask(ji,j2+1,jk) 
     838               END DO  
     839            END DO 
     840         ENDIF 
     841         ! 
     842      ENDIF 
     843      !  
     844   END SUBROUTINE correct_v_bdy 
    290845 
    291846 
    … …  
    298853      LOGICAL                         , INTENT(in   ) ::   before 
    299854      !!  
    300       INTEGER ::   ji, jj, jk 
    301       REAL(wp)::   zrhoy, zcorr 
     855      INTEGER  :: ji, jj, jk 
     856      REAL(wp) :: zrhoy 
     857      REAL(wp) :: zcorr 
    302858      !!--------------------------------------------- 
    303859      ! 
    … …  
    312868         DO jj=j1,j2 
    313869            DO ji=i1,i2 
    314                tabres(ji,jj) =  tabres(ji,jj) * r1_hu_n(ji,jj) * r1_e2u(ji,jj)   
     870               tabres(ji,jj) =  tabres(ji,jj) * r1_e2u(ji,jj)   
    315871               !     
    316872               ! Update "now" 3d velocities: 
    … …  
    319875                  spgu(ji,jj) = spgu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) 
    320876               END DO 
    321                spgu(ji,jj) = spgu(ji,jj) * r1_hu_n(ji,jj) 
    322877               ! 
    323                zcorr = tabres(ji,jj) - spgu(ji,jj) 
     878               zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu_n(ji,jj) 
    324879               DO jk=1,jpkm1               
    325880                  un(ji,jj,jk) = un(ji,jj,jk) + zcorr * umask(ji,jj,jk)            
    … …  
    329884               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 
    330885                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    331                      zcorr = tabres(ji,jj) - un_b(ji,jj) 
     886                     zcorr = (tabres(ji,jj) - un_b(ji,jj) * hu_a(ji,jj)) * r1_hu_b(ji,jj) 
    332887                     ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1) 
    333888                  END IF 
    334                ENDIF              
    335                un_b(ji,jj) = tabres(ji,jj) * umask(ji,jj,1) 
     889               ENDIF     
     890               un_b(ji,jj) = tabres(ji,jj) * r1_hu_n(ji,jj) * umask(ji,jj,1) 
    336891               !        
    337892               ! Correct "before" velocities to hold correct bt component: 
    … …  
    340895                  spgu(ji,jj) = spgu(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) 
    341896               END DO 
    342                spgu(ji,jj) = spgu(ji,jj) * r1_hu_b(ji,jj) 
    343897               ! 
    344                zcorr = ub_b(ji,jj) - spgu(ji,jj) 
     898               zcorr = ub_b(ji,jj) - spgu(ji,jj) * r1_hu_b(ji,jj) 
    345899               DO jk=1,jpkm1               
    346900                  ub(ji,jj,jk) = ub(ji,jj,jk) + zcorr * umask(ji,jj,jk)            
    … …  
    349903            END DO 
    350904         END DO 
     905         ! 
     906         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     907            ub_b(i1:i2,j1:j2)  = un_b(i1:i2,j1:j2) 
     908         ENDIF 
    351909      ENDIF 
    352910      ! 
    … …  
    362920      LOGICAL                         , INTENT(in   ) ::   before 
    363921      !  
    364       INTEGER :: ji, jj, jk 
     922      INTEGER  :: ji, jj, jk 
    365923      REAL(wp) :: zrhox, zcorr 
    366924      !!---------------------------------------------------------------------- 
    … …  
    376934         DO jj=j1,j2 
    377935            DO ji=i1,i2 
    378                tabres(ji,jj) =  tabres(ji,jj) * r1_hv_n(ji,jj) * r1_e1v(ji,jj)   
     936               tabres(ji,jj) =  tabres(ji,jj) * r1_e1v(ji,jj)   
    379937               !     
    380938               ! Update "now" 3d velocities: 
    … …  
    383941                  spgv(ji,jj) = spgv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) 
    384942               END DO 
    385                spgv(ji,jj) = spgv(ji,jj) * r1_hv_n(ji,jj) 
    386943               ! 
    387                zcorr = tabres(ji,jj) - spgv(ji,jj) 
     944               zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv_n(ji,jj) 
    388945               DO jk=1,jpkm1               
    389946                  vn(ji,jj,jk) = vn(ji,jj,jk) + zcorr * vmask(ji,jj,jk)            
    … …  
    393950               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 
    394951                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    395                      zcorr = tabres(ji,jj) - vn_b(ji,jj) 
     952                     zcorr = (tabres(ji,jj) - vn_b(ji,jj) * hv_a(ji,jj)) * r1_hv_b(ji,jj) 
    396953                     vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1) 
    397954                  END IF 
    398955               ENDIF               
    399                vn_b(ji,jj) = tabres(ji,jj) * vmask(ji,jj,1) 
     956               vn_b(ji,jj) = tabres(ji,jj) * r1_hv_n(ji,jj) * vmask(ji,jj,1) 
    400957               !        
    401958               ! Correct "before" velocities to hold correct bt component: 
    … …  
    404961                  spgv(ji,jj) = spgv(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) 
    405962               END DO 
    406                spgv(ji,jj) = spgv(ji,jj) * r1_hv_b(ji,jj) 
    407963               ! 
    408                zcorr = vb_b(ji,jj) - spgv(ji,jj) 
     964               zcorr = vb_b(ji,jj) - spgv(ji,jj) * r1_hv_b(ji,jj) 
    409965               DO jk=1,jpkm1               
    410966                  vb(ji,jj,jk) = vb(ji,jj,jk) + zcorr * vmask(ji,jj,jk)            
    … …  
    413969            END DO 
    414970         END DO 
     971         ! 
     972         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     973            vb_b(i1:i2,j1:j2)  = vn_b(i1:i2,j1:j2) 
     974         ENDIF 
     975         ! 
    415976      ENDIF 
    416977      !  
    … …  
    436997         END DO 
    437998      ELSE 
    438          IF( .NOT.ln_dynspg_ts .OR. ( ln_dynspg_ts .AND. .NOT.ln_bt_fw ) ) THEN 
    439             IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
    440                DO jj=j1,j2 
    441                   DO ji=i1,i2 
    442                      sshb(ji,jj) =   sshb(ji,jj) & 
    443                            & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1) 
    444                   END DO 
    445                END DO 
    446             ENDIF 
     999         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
     1000            DO jj=j1,j2 
     1001               DO ji=i1,i2 
     1002                  sshb(ji,jj) =   sshb(ji,jj) & 
     1003                        & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1) 
     1004               END DO 
     1005            END DO 
    4471006         ENDIF 
    4481007         ! 
    … …  
    4521011            END DO 
    4531012         END DO 
     1013         ! 
     1014         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     1015            sshb(i1:i2,j1:j2)  = sshn(i1:i2,j1:j2) 
     1016         ENDIF 
     1017         ! 
     1018 
    4541019      ENDIF 
    4551020      ! 
    … …  
    4651030      LOGICAL                            , INTENT(in) ::   before 
    4661031      !! 
    467       INTEGER ::   ji, jj 
    468       REAL(wp)::   zrhoy 
    469       !!---------------------------------------------------------------------- 
     1032      INTEGER :: ji, jj 
     1033      REAL(wp) :: zrhoy, za1, zcor 
     1034      !!--------------------------------------------- 
    4701035      ! 
    4711036      IF (before) THEN 
    … …  
    4781043         tabres = zrhoy * tabres 
    4791044      ELSE 
     1045         !  
     1046         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2) 
     1047         ! 
     1048         za1 = 1._wp / REAL(Agrif_rhot(), wp) 
    4801049         DO jj=j1,j2 
    4811050            DO ji=i1,i2 
    482                ub2_b(ji,jj) = tabres(ji,jj) / e2u(ji,jj) 
     1051               zcor=tabres(ji,jj) - ub2_b(ji,jj) 
     1052               ! Update time integrated fluxes also in case of multiply nested grids: 
     1053               ub2_i_b(ji,jj) = ub2_i_b(ji,jj) + za1 * zcor  
     1054               ! Update corrective fluxes: 
     1055               un_bf(ji,jj)  = un_bf(ji,jj) + zcor 
     1056               ! Update half step back fluxes: 
     1057               ub2_b(ji,jj) = tabres(ji,jj) 
    4831058            END DO 
    4841059         END DO 
    … …  
    4871062   END SUBROUTINE updateub2b 
    4881063 
     1064   SUBROUTINE reflux_sshu( tabres, i1, i2, j1, j2, before, nb, ndir ) 
     1065      !!--------------------------------------------- 
     1066      !!          *** ROUTINE reflux_sshu *** 
     1067      !!--------------------------------------------- 
     1068      INTEGER, INTENT(in) :: i1, i2, j1, j2 
     1069      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres 
     1070      LOGICAL, INTENT(in) :: before 
     1071      INTEGER, INTENT(in) :: nb, ndir 
     1072      !! 
     1073      LOGICAL :: western_side, eastern_side  
     1074      INTEGER :: ji, jj 
     1075      REAL(wp) :: zrhoy, za1, zcor 
     1076      !!--------------------------------------------- 
     1077      ! 
     1078      IF (before) THEN 
     1079         zrhoy = Agrif_Rhoy() 
     1080         DO jj=j1,j2 
     1081            DO ji=i1,i2 
     1082               tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj) 
     1083            END DO 
     1084         END DO 
     1085         tabres = zrhoy * tabres 
     1086      ELSE 
     1087         !  
     1088         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2) 
     1089         ! 
     1090         western_side  = (nb == 1).AND.(ndir == 1) 
     1091         eastern_side  = (nb == 1).AND.(ndir == 2) 
     1092         ! 
     1093         IF (western_side) THEN 
     1094            DO jj=j1,j2 
     1095               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 
     1098            END DO 
     1099         ENDIF 
     1100         IF (eastern_side) THEN 
     1101            DO jj=j1,j2 
     1102               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 
     1105            END DO 
     1106         ENDIF 
     1107         ! 
     1108      ENDIF 
     1109      ! 
     1110   END SUBROUTINE reflux_sshu 
    4891111 
    4901112   SUBROUTINE updatevb2b( tabres, i1, i2, j1, j2, before ) 
    … …  
    4961118      LOGICAL                         , INTENT(in   ) ::   before 
    4971119      !! 
    498       INTEGER ::   ji, jj 
    499       REAL(wp)::   zrhox 
    500       !!---------------------------------------------------------------------- 
     1120      INTEGER :: ji, jj 
     1121      REAL(wp) :: zrhox, za1, zcor 
     1122      !!--------------------------------------------- 
    5011123      ! 
    5021124      IF( before ) THEN 
    … …  
    5091131         tabres = zrhox * tabres 
    5101132      ELSE 
     1133         !  
     1134         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2) 
     1135         ! 
     1136         za1 = 1._wp / REAL(Agrif_rhot(), wp) 
    5111137         DO jj=j1,j2 
    5121138            DO ji=i1,i2 
    513                vb2_b(ji,jj) = tabres(ji,jj) / e1v(ji,jj) 
     1139               zcor=tabres(ji,jj) - vb2_b(ji,jj) 
     1140               ! Update time integrated fluxes also in case of multiply nested grids: 
     1141               vb2_i_b(ji,jj) = vb2_i_b(ji,jj) + za1 * zcor  
     1142               ! Update corrective fluxes: 
     1143               vn_bf(ji,jj)  = vn_bf(ji,jj) + zcor 
     1144               ! Update half step back fluxes: 
     1145               vb2_b(ji,jj) = tabres(ji,jj) 
    5141146            END DO 
    5151147         END DO 
    … …  
    5181150   END SUBROUTINE updatevb2b 
    5191151 
     1152   SUBROUTINE reflux_sshv( tabres, i1, i2, j1, j2, before, nb, ndir ) 
     1153      !!--------------------------------------------- 
     1154      !!          *** ROUTINE reflux_sshv *** 
     1155      !!--------------------------------------------- 
     1156      INTEGER, INTENT(in) :: i1, i2, j1, j2 
     1157      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres 
     1158      LOGICAL, INTENT(in) :: before 
     1159      INTEGER, INTENT(in) :: nb, ndir 
     1160      !! 
     1161      LOGICAL :: southern_side, northern_side  
     1162      INTEGER :: ji, jj 
     1163      REAL(wp) :: zrhox, za1, zcor 
     1164      !!--------------------------------------------- 
     1165      ! 
     1166      IF (before) THEN 
     1167         zrhox = Agrif_Rhox() 
     1168         DO jj=j1,j2 
     1169            DO ji=i1,i2 
     1170               tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj)  
     1171            END DO 
     1172         END DO 
     1173         tabres = zrhox * tabres 
     1174      ELSE 
     1175         !  
     1176         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2) 
     1177         ! 
     1178         southern_side = (nb == 2).AND.(ndir == 1) 
     1179         northern_side = (nb == 2).AND.(ndir == 2) 
     1180         ! 
     1181         IF (southern_side) THEN 
     1182            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 
     1186            END DO 
     1187         ENDIF 
     1188         IF (northern_side) THEN                
     1189            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 
     1193            END DO 
     1194         ENDIF 
     1195         !  
     1196      ENDIF 
     1197      ! 
     1198   END SUBROUTINE reflux_sshv 
    5201199 
    5211200   SUBROUTINE update_scales( tabres, i1, i2, j1, j2, k1, k2, n1,n2, before ) 
    … …  
    6191298   END SUBROUTINE updateAVM 
    6201299 
     1300   SUBROUTINE updatee3t(ptab_dum, i1, i2, j1, j2, k1, k2, before ) 
     1301      !!--------------------------------------------- 
     1302      !!           *** ROUTINE updatee3t *** 
     1303      !!--------------------------------------------- 
     1304      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ptab_dum 
     1305      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2 
     1306      LOGICAL, INTENT(in) :: before 
     1307      ! 
     1308      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptab 
     1309      INTEGER :: ji,jj,jk 
     1310      REAL(wp) :: zcoef 
     1311      !!--------------------------------------------- 
     1312      ! 
     1313      IF (.NOT.before) THEN 
     1314         ! 
     1315         ALLOCATE(ptab(i1:i2,j1:j2,1:jpk))  
     1316         ! 
     1317         ! Update e3t from ssh (z* case only) 
     1318         DO jk = 1, jpkm1 
     1319            DO jj=j1,j2 
     1320               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))) 
     1323               END DO 
     1324            END DO 
     1325         END DO 
     1326         ! 
     1327         ! 1) Updates at BEFORE time step: 
     1328         ! ------------------------------- 
     1329         ! 
     1330         ! Save "old" scale factor (prior update) for subsequent asselin correction 
     1331         ! 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) 
     1336 
     1337         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN 
     1338            DO jk = 1, jpkm1 
     1339               DO jj=j1,j2 
     1340                  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) ) 
     1343                  END DO 
     1344               END DO 
     1345            END DO 
     1346            ! 
     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) 
     1350            ! 
     1351            DO jk = 2, jpk 
     1352               DO jj = j1,j2 
     1353                  DO ji = i1,i2             
     1354                     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))  
     1362                  END DO 
     1363               END DO 
     1364            END DO 
     1365            ! 
     1366         ENDIF         
     1367         ! 
     1368         ! 2) Updates at NOW time step: 
     1369         ! ---------------------------- 
     1370         ! 
     1371         ! Update vertical scale factor at T-points: 
     1372         e3t_n(i1:i2,j1:j2,1:jpkm1) = ptab(i1:i2,j1:j2,1:jpkm1) 
     1373         ! 
     1374         ! Update total depth: 
     1375         ht_n(i1:i2,j1:j2) = 0._wp 
     1376         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) 
     1378         END DO 
     1379         ! 
     1380         ! 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 
     1385         ! 
     1386         DO jk = 2, jpk 
     1387            DO jj = j1,j2 
     1388               DO ji = i1,i2             
     1389               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 
     1396               END DO 
     1397            END DO 
     1398         END DO 
     1399         ! 
     1400         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) 
     1405         ENDIF 
     1406         ! 
     1407         DEALLOCATE(ptab) 
     1408      ENDIF 
     1409      ! 
     1410   END SUBROUTINE updatee3t 
     1411 
    6211412#else 
    6221413   !!---------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.