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 14086 – NEMO

Changeset 14086


Ignore:
Timestamp:
2020-12-04T12:37:14+01:00 (4 years ago)
Author:
cetlod
Message:

Adding AGRIF branches into the trunk

Location:
NEMO/trunk
Files:
53 edited
29 copied

Legend:

Unmodified
Added
Removed

  • NEMO/trunk/cfgs/AGRIF_DEMO/EXPREF/1_context_nemo.xml

    r13476 r14086  
    44==============================================================================================  
    55--> 
    6 <context id="1_nemo"> 
     6<context id="nemo"> 
    77    <!-- $id$ --> 
    88    <variable_definition> 
    … …  
    2222    <field_definition src="./field_def_nemo-oce.xml"/>   <!--  NEMO ocean dynamics                     --> 
    2323    <field_definition src="./field_def_nemo-ice.xml"/>   <!--  NEMO ocean sea ice                      --> 
     24    <field_definition src="./field_def_nemo-pisces.xml"/>   <!--  NEMO ocean biogeochemical                      --> 
     25    <field_definition src="./field_def_nemo-innerttrc.xml"/> <!--  NEMO ocean inert passive tracer           --> 
    2426 
    2527 
    … …  
    2729    <file_definition src="./file_def_nemo-oce.xml"/>     <!--  NEMO ocean dynamics                     --> 
    2830    <file_definition src="./file_def_nemo-ice.xml"/>     <!--  NEMO ocean sea ice                      --> 
     31    <file_definition src="./file_def_nemo-innerttrc.xml"/> <!--  NEMO ocean inert passive tracer           --> 
    2932 
    3033<!-- Axis definition --> 

  • NEMO/trunk/cfgs/AGRIF_DEMO/EXPREF/2_context_nemo.xml

    r13476 r14086  
    44==============================================================================================  
    55--> 
    6 <context id="2_nemo"> 
     6<context id="nemo"> 
    77    <!-- $id$ --> 
    88    <variable_definition> 
    … …  
    2222    <field_definition src="./field_def_nemo-oce.xml"/>   <!--  NEMO ocean dynamics                     --> 
    2323    <field_definition src="./field_def_nemo-ice.xml"/>   <!--  NEMO ocean sea ice                      --> 
     24    <field_definition src="./field_def_nemo-pisces.xml"/>   <!--  NEMO ocean biogeochemical                      --> 
     25    <field_definition src="./field_def_nemo-innerttrc.xml"/> <!--  NEMO ocean inert passive tracer           --> 
    2426 
    2527 
    … …  
    2729    <file_definition src="./file_def_nemo-oce.xml"/>     <!--  NEMO ocean dynamics                     --> 
    2830    <file_definition src="./file_def_nemo-ice.xml"/>     <!--  NEMO ocean sea ice                      --> 
     31    <file_definition src="./file_def_nemo-innerttrc.xml"/> <!--  NEMO ocean inert passive tracer           --> 
    2932 
    3033<!-- Axis definition --> 

  • NEMO/trunk/cfgs/AGRIF_DEMO/EXPREF/2_namelist_cfg

    r13558 r14086  
    158158!----------------------------------------------------------------------- 
    159159   ln_spc_dyn    = .true.  !  use 0 as special value for dynamics 
    160    ln_chk_bathy  = .true. !  =T  check the parent bathymetry 
     160   ln_chk_bathy  = .false. !  =T  check the parent bathymetry 
    161161/ 
    162162!----------------------------------------------------------------------- 

  • NEMO/trunk/cfgs/AGRIF_DEMO/EXPREF/3_context_nemo.xml

    r13476 r14086  
    44==============================================================================================  
    55--> 
    6 <context id="3_nemo"> 
     6<context id="nemo"> 
    77    <!-- $id$ --> 
    88    <variable_definition> 
    … …  
    2222    <field_definition src="./field_def_nemo-oce.xml"/>   <!--  NEMO ocean dynamics                     --> 
    2323    <field_definition src="./field_def_nemo-ice.xml"/>   <!--  NEMO ocean sea ice                      --> 
     24    <field_definition src="./field_def_nemo-pisces.xml"/>   <!--  NEMO ocean biogeochemical                      --> 
     25    <field_definition src="./field_def_nemo-innerttrc.xml"/> <!--  NEMO ocean inert passive tracer           --> 
    2426 
    2527 
    … …  
    2729    <file_definition src="./file_def_nemo-oce.xml"/>     <!--  NEMO ocean dynamics                     --> 
    2830    <file_definition src="./file_def_nemo-ice.xml"/>     <!--  NEMO ocean sea ice                      --> 
     31    <file_definition src="./file_def_nemo-innerttrc.xml"/> <!--  NEMO ocean inert passive tracer           --> 
    2932 
    3033<!-- Axis definition --> 

  • NEMO/trunk/cfgs/AGRIF_DEMO/EXPREF/3_namelist_cfg

    r13558 r14086  
    158158!----------------------------------------------------------------------- 
    159159   ln_spc_dyn    = .true.  !  use 0 as special value for dynamics 
    160    ln_chk_bathy  = .true.  !  =T  check the parent bathymetry 
     160   ln_chk_bathy  = .false.  !  =T  check the parent bathymetry 
    161161/ 
    162162!----------------------------------------------------------------------- 

  • NEMO/trunk/cfgs/AGRIF_DEMO/EXPREF/context_nemo.xml

    r13476 r14086  
    2222    <field_definition src="./field_def_nemo-oce.xml"/>   <!--  NEMO ocean dynamics                     --> 
    2323    <field_definition src="./field_def_nemo-ice.xml"/>   <!--  NEMO ocean sea ice                      --> 
     24    <field_definition src="./field_def_nemo-pisces.xml"/>   <!--  NEMO ocean biogeochemical                      --> 
     25    <field_definition src="./field_def_nemo-innerttrc.xml"/> <!--  NEMO ocean inert passive tracer           --> 
    2426 
    2527 
    … …  
    2729    <file_definition src="./file_def_nemo-oce.xml"/>     <!--  NEMO ocean dynamics                     --> 
    2830    <file_definition src="./file_def_nemo-ice.xml"/>     <!--  NEMO ocean sea ice                      --> 
     31    <file_definition src="./file_def_nemo-innerttrc.xml"/> <!--  NEMO ocean inert passive tracer           --> 
    2932 
    3033<!-- Axis definition --> 

  • NEMO/trunk/cfgs/AGRIF_DEMO/cpp_AGRIF_DEMO.fcm

    r12208 r14086  
    1 bld::tool::fppkeys   key_si3 key_iomput key_mpp_mpi key_agrif 
     1bld::tool::fppkeys   key_si3 key_top key_iomput key_mpp_mpi key_agrif 

  • NEMO/trunk/cfgs/ref_cfgs.txt

    r14072 r14086  
    1 AGRIF_DEMO OCE ICE NST 
     1AGRIF_DEMO OCE ICE TOP NST 
    22AMM12 OCE 
    33C1D_PAPA OCE 

  • NEMO/trunk/src/ICE/iceistate.F90

    r14072 r14086  
    1818   USE oce            ! dynamics and tracers variables 
    1919   USE dom_oce        ! ocean domain 
    20    USE sbc_oce , ONLY : sst_m, sss_m, ln_ice_embd 
     20   USE sbc_oce , ONLY : sst_m, sss_m, ln_ice_embd  
    2121   USE sbc_ice , ONLY : tn_ice, snwice_mass, snwice_mass_b 
    2222   USE eosbn2         ! equation of state 
    … …  
    3939# if defined key_agrif 
    4040   USE agrif_oce 
    41    USE agrif_ice 
    42    USE agrif_ice_interp 
    43 # endif 
     41   USE agrif_ice_interp  
     42# endif    
    4443 
    4544   IMPLICIT NONE 
    … …  
    9190      !! 
    9291      !! ** Method  :   This routine will put some ice where ocean 
    93       !!                is at the freezing point, then fill in ice 
    94       !!                state variables using prescribed initial 
    95       !!                values in the namelist 
     92      !!                is at the freezing point, then fill in ice  
     93      !!                state variables using prescribed initial  
     94      !!                values in the namelist             
    9695      !! 
    9796      !! ** Steps   :   1) Set initial surface and basal temperatures 
    … …  
    103102      !!              where there is no ice 
    104103      !!-------------------------------------------------------------------- 
    105       INTEGER, INTENT(in) :: kt            ! time step 
     104      INTEGER, INTENT(in) :: kt            ! time step  
    106105      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices 
    107106      ! 
    … …  
    129128      ! basal temperature (considered at freezing point)   [Kelvin] 
    130129      CALL eos_fzp( sss_m(:,:), t_bo(:,:) ) 
    131       t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 
     130      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1)  
    132131      ! 
    133132      ! surface temperature and conductivity 
    … …  
    154153      e_i (:,:,:,:) = 0._wp 
    155154      e_s (:,:,:,:) = 0._wp 
    156  
     155       
    157156      ! general fields 
    158157      a_i (:,:,:) = 0._wp 
    … …  
    184183      IF( ln_iceini ) THEN 
    185184         ! 
    186          IF( Agrif_Root() ) THEN 
     185#if defined key_agrif 
     186         IF ( ( Agrif_Root() ).OR.(.NOT.ln_init_chfrpar ) ) THEN 
     187#endif 
    187188            !                             !---------------! 
    188189            IF( nn_iceini_file == 1 )THEN ! Read a file   ! 
    … …  
    229230               IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) & 
    230231                  &     si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) & ! rn_apd = pond fraction => rn_apnd * a_i = pond conc. 
    231                   &                              * si(jp_ati)%fnow(:,:,1) 
     232                  &                              * si(jp_ati)%fnow(:,:,1)  
    232233               ! 
    233234               ! pond depth 
    … …  
    248249               ! 
    249250               ! change the switch for the following 
    250                WHERE( zat_i_ini(:,:) > 0._wp )   ;   zswitch(:,:) = tmask(:,:,1) 
     251               WHERE( zat_i_ini(:,:) > 0._wp )   ;   zswitch(:,:) = tmask(:,:,1)  
    251252               ELSEWHERE                         ;   zswitch(:,:) = 0._wp 
    252253               END WHERE 
    … …  
    256257               !                          !---------------! 
    257258               ! no ice if (sst - Tfreez) >= thresold 
    258                WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst )   ;   zswitch(:,:) = 0._wp 
     259               WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst )   ;   zswitch(:,:) = 0._wp  
    259260               ELSEWHERE                                                                    ;   zswitch(:,:) = tmask(:,:,1) 
    260261               END WHERE 
    … …  
    269270                  zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:) 
    270271                  ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:) 
    271                   zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 
     272                  zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.  
    272273                  zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:) 
    273274                  zhlid_ini(:,:) = rn_hld_ini_n * zswitch(:,:) 
    … …  
    295296               zhlid_ini(:,:) = 0._wp 
    296297            ENDIF 
    297  
     298             
    298299            IF ( .NOT.ln_pnd_lids ) THEN 
    299300               zhlid_ini(:,:) = 0._wp 
    300301            ENDIF 
    301  
     302             
    302303            !----------------! 
    303304            ! 3) fill fields ! 
    … …  
    323324            CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti)  , zhpnd_ini ) 
    324325            CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d(1:npti)  , zhlid_ini ) 
    325  
     326             
    326327            ! allocate temporary arrays 
    327328            ALLOCATE( zhi_2d (npti,jpl), zhs_2d (npti,jpl), zai_2d (npti,jpl), & 
    … …  
    377378            DO jl = 1, jpl 
    378379               DO_3D( 1, 1, 1, 1, 1, nlay_i ) 
    379                   t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl) 
     380                  t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl)  
    380381                  ztmelts          = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K 
    381382                  e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * & 
    … …  
    385386               END_3D 
    386387            END DO 
    387  
    388 #if  defined key_agrif 
     388            ! 
     389#if defined key_agrif 
    389390         ELSE 
    390  
    391             Agrif_SpecialValue    = -9999. 
    392             Agrif_UseSpecialValue = .TRUE. 
    393             CALL Agrif_init_variable(tra_iceini_id,procname=interp_tra_ice) 
    394             use_sign_north = .TRUE. 
    395             sign_north = -1. 
    396             CALL Agrif_init_variable(u_iceini_id  ,procname=interp_u_ice) 
    397             CALL Agrif_init_variable(v_iceini_id  ,procname=interp_v_ice) 
    398             Agrif_SpecialValue    = 0._wp 
    399             use_sign_north = .FALSE. 
    400             Agrif_UseSpecialValue = .FALSE. 
    401         ! lbc ???? 
    402    ! Here we know : a_i, v_i, v_s, sv_i, oa_i, a_ip, v_ip, v_il, t_su, e_s, e_i 
    403             CALL ice_var_glo2eqv 
    404             CALL ice_var_zapsmall 
    405             CALL ice_var_agg(2) 
     391            CALL  agrif_istate_ice 
     392         ENDIF 
    406393#endif 
    407          ENDIF ! Agrif_Root 
    408          ! 
    409394         ! Melt ponds 
    410395         WHERE( a_i > epsi10 )   ;   a_ip_eff(:,:,:) = a_ip(:,:,:) / a_i(:,:,:) 
    … …  
    413398         v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:) 
    414399         v_il(:,:,:) = h_il(:,:,:) * a_ip(:,:,:) 
    415  
     400          
    416401         ! specific temperatures for coupled runs 
    417402         tn_ice(:,:,:) = t_su(:,:,:) 
    … …  
    423408            WHERE( at_i(:,:) > rn_amax_2d(:,:) )   a_i(:,:,jl) = a_i(:,:,jl) * rn_amax_2d(:,:) / at_i(:,:) 
    424409         END DO 
    425          at_i(:,:) = SUM( a_i, dim=3 ) 
     410        at_i(:,:) = SUM( a_i, dim=3 ) 
    426411         ! 
    427412      ENDIF ! ln_iceini 
    … …  
    456441      !!------------------------------------------------------------------- 
    457442      !!                   ***  ROUTINE ice_istate_init  *** 
    458       !! 
    459       !! ** Purpose :   Definition of initial state of the ice 
    460       !! 
    461       !! ** Method  :   Read the namini namelist and check the parameter 
     443      !!         
     444      !! ** Purpose :   Definition of initial state of the ice  
     445      !! 
     446      !! ** Method  :   Read the namini namelist and check the parameter  
    462447      !!              values called at the first timestep (nit000) 
    463448      !! 
    … …  
    500485         WRITE(numout,*) '      max ocean temp. above Tfreeze with initial ice   rn_thres_sst   = ', rn_thres_sst 
    501486         IF( ln_iceini .AND. nn_iceini_file == 0 ) THEN 
    502             WRITE(numout,*) '      initial snw thickness in the north-south         rn_hts_ini     = ', rn_hts_ini_n,rn_hts_ini_s 
     487            WRITE(numout,*) '      initial snw thickness in the north-south         rn_hts_ini     = ', rn_hts_ini_n,rn_hts_ini_s  
    503488            WRITE(numout,*) '      initial ice thickness in the north-south         rn_hti_ini     = ', rn_hti_ini_n,rn_hti_ini_s 
    504489            WRITE(numout,*) '      initial ice concentr  in the north-south         rn_ati_ini     = ', rn_ati_ini_n,rn_ati_ini_s 

  • NEMO/trunk/src/NST/agrif_ice_interp.F90

    r13472 r14086  
    2525   USE agrif_oce 
    2626   USE phycst , ONLY: rt0 
     27   USE icevar 
     28   USE sbc_ice, ONLY : tn_ice 
    2729    
    2830   IMPLICIT NONE 
    … …  
    3032 
    3133   PUBLIC   agrif_interp_ice   ! called by agrif_user.F90 
    32    PUBLIC   interp_tra_ice, interp_u_ice, interp_v_ice  ! called by iceistate.F90 
     34   PUBLIC   agrif_istate_ice   ! called by icerst.F90 
    3335 
    3436   !!---------------------------------------------------------------------- 
    … …  
    3941 
    4042CONTAINS 
     43 
     44   SUBROUTINE agrif_istate_ice 
     45      !!----------------------------------------------------------------------- 
     46      !!                 *** ROUTINE agrif_istate_ice  *** 
     47      !! 
     48      !!  ** Method  : Set initial ice fields from parent grid 
     49      !! 
     50      !!----------------------------------------------------------------------- 
     51      IF(lwp) WRITE(numout,*) ' ' 
     52      IF(lwp) WRITE(numout,*) 'Agrif_istate_ice : interp child ice initial state from parent' 
     53      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~' 
     54      IF(lwp) WRITE(numout,*) ' ' 
     55 
     56      ! Set a_i, v_i, v_s, sv_i, oa_i, a_ip, v_ip, t_su, e_s, e_i: 
     57      Agrif_SpecialValue    = -9999. 
     58      Agrif_UseSpecialValue = .TRUE. 
     59      CALL Agrif_init_variable(tra_iceini_id,procname=interp_tra_ice) 
     60      ! 
     61      ! Set u_ice, v_ice: 
     62      use_sign_north = .TRUE. 
     63      sign_north = -1. 
     64      Agrif_UseSpecialValue = .TRUE. 
     65      CALL Agrif_init_variable(u_iceini_id  ,procname=interp_u_ice) 
     66      CALL Agrif_init_variable(v_iceini_id  ,procname=interp_v_ice) 
     67      Agrif_SpecialValue = 0._wp 
     68      use_sign_north = .FALSE. 
     69      Agrif_UseSpecialValue = .FALSE. 
     70      ! lbc ???? 
     71      ! JC: do we really need the 3 lines below ? 
     72      CALL ice_var_glo2eqv 
     73      CALL ice_var_zapsmall 
     74      CALL ice_var_agg(2) 
     75 
     76      ! Melt ponds 
     77      WHERE( a_i > epsi10 ) 
     78         a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:) 
     79      ELSEWHERE 
     80         a_ip_frac(:,:,:) = 0._wp 
     81      END WHERE 
     82      WHERE( a_ip > 0._wp )       ! ???????     
     83         h_ip(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:) 
     84      ELSEWHERE 
     85         h_ip(:,:,:) = 0._wp 
     86      END WHERE    
     87 
     88      tn_ice(:,:,:) = t_su(:,:,:) 
     89      t1_ice(:,:,:) = t_i (:,:,1,:)  
     90 
     91   END SUBROUTINE agrif_istate_ice 
    4192 
    4293   SUBROUTINE agrif_interp_ice( cd_type, kiter, kitermax ) 

  • NEMO/trunk/src/NST/agrif_oce.F90

    r13286 r14086  
    2323   LOGICAL , PUBLIC ::   ln_spc_dyn    = .FALSE.   !: use zeros (.false.) or not (.true.) in 
    2424                                                   !: bdys dynamical fields interpolation 
    25    REAL(wp), PUBLIC ::   rn_sponge_tra = 2800.     !: sponge coeff. for tracers 
    26    REAL(wp), PUBLIC ::   rn_sponge_dyn = 2800.     !: sponge coeff. for dynamics 
     25   LOGICAL , PUBLIC ::   ln_vert_remap = .FALSE.   !: use vertical remapping 
     26   REAL(wp), PUBLIC ::   rn_sponge_tra = 0.002     !: sponge coeff. for tracers 
     27   REAL(wp), PUBLIC ::   rn_sponge_dyn = 0.002     !: sponge coeff. for dynamics 
    2728   REAL(wp), PUBLIC ::   rn_trelax_tra = 0.01      !: time relaxation parameter for tracers 
    2829   REAL(wp), PUBLIC ::   rn_trelax_dyn = 0.01      !: time relaxation parameter for momentum 
    … …  
    3031   ! 
    3132   INTEGER , PUBLIC, PARAMETER ::   nn_sponge_len = 2  !: Sponge width (in number of parent grid points) 
     33   INTEGER , PUBLIC, PARAMETER ::   nn_shift_bar = 0   !: nb of coarse grid points by which we shift 2d interface 
    3234 
    3335   LOGICAL , PUBLIC :: spongedoneT = .FALSE.       !: tracer   sponge layer indicator 
    … …  
    3537   LOGICAL , PUBLIC :: lk_agrif_fstep = .TRUE.     !: if true: first step 
    3638   LOGICAL , PUBLIC :: lk_agrif_debug = .FALSE.    !: if true: print debugging info 
    37  
     39   LOGICAL , PUBLIC :: lk_tint2d_notinterp = .FALSE. !: if true, no time interp 
    3840   LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_tsn 
    3941# if defined key_top 
    … …  
    4648   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fspu, fspv !: sponge arrays 
    4749   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fspt, fspf !:   "      " 
     50   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fspu_2d,fspv_2d  !: sponge arrays (2d mode) 
     51   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fspt_2d, fspf_2d !:   "       "     "   " 
    4852 
    4953   ! Barotropic arrays used to store open boundary data during time-splitting loop: 
    … …  
    5155   INTEGER , PUBLIC,              SAVE                 ::  Kbb_a, Kmm_a, Krhs_a   !: AGRIF module-specific copies of time-level indices 
    5256 
    53    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht0_parent, hu0_parent, hv0_parent 
    54    INTEGER,  PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt_parent, mbku_parent, mbkv_parent 
     57   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: ht0_parent, hu0_parent, hv0_parent 
     58   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t0_parent, e3u0_parent, e3v0_parent 
     59   INTEGER,  PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: mbkt_parent, mbku_parent, mbkv_parent 
    5560 
    56    INTEGER, PUBLIC :: tsn_id                                                  ! AGRIF profile for tracers interpolation and update 
     61 
     62   INTEGER, PUBLIC :: ts_interp_id, ts_update_id                              ! AGRIF profile for tracers interpolation and update 
    5763   INTEGER, PUBLIC :: un_interp_id, vn_interp_id                              ! AGRIF profiles for interpolations 
    5864   INTEGER, PUBLIC :: un_update_id, vn_update_id                              ! AGRIF profiles for udpates 
    59    INTEGER, PUBLIC :: tsn_sponge_id, un_sponge_id, vn_sponge_id               ! AGRIF profiles for sponge layers 
     65   INTEGER, PUBLIC :: ts_sponge_id, un_sponge_id, vn_sponge_id                ! AGRIF profiles for sponge layers (3d) 
     66   INTEGER, PUBLIC :: unb_sponge_id, vnb_sponge_id                            ! AGRIF profiles for sponge layers (2d) 
    6067   INTEGER, PUBLIC :: tsini_id, uini_id, vini_id, sshini_id                   ! AGRIF profile for initialization 
    6168# if defined key_top 
    6269   INTEGER, PUBLIC :: trn_id, trn_sponge_id 
    6370# endif   
    64    INTEGER, PUBLIC :: unb_id, vnb_id, ub2b_interp_id, vb2b_interp_id 
    65    INTEGER, PUBLIC :: ub2b_update_id, vb2b_update_id 
    66    INTEGER, PUBLIC :: e3t_id, e1u_id, e2v_id, sshn_id 
     71   INTEGER, PUBLIC :: unb_interp_id, vnb_interp_id, ub2b_interp_id, vb2b_interp_id 
     72   INTEGER, PUBLIC :: ub2b_update_id, vb2b_update_id, unb_update_id, vnb_update_id 
     73   INTEGER, PUBLIC :: ub2b_cor_id, vb2b_cor_id 
     74   INTEGER, PUBLIC :: e3t_id, sshn_id 
    6775   INTEGER, PUBLIC :: scales_t_id 
    6876   INTEGER, PUBLIC :: avt_id, avm_id, en_id                ! TKE related identificators 
    69    INTEGER, PUBLIC :: mbkt_id, ht0_id 
     77   INTEGER, PUBLIC :: mbkt_id, ht0_id, e3t0_interp_id 
    7078   INTEGER, PUBLIC :: glamt_id, gphit_id 
     79   INTEGER, PUBLIC :: batupd_id 
    7180   INTEGER, PUBLIC :: kindic_agr 
    7281 
    … …  
    7483!$AGRIF_DO_NOT_TREAT 
    7584   LOGICAL, PUBLIC :: use_sign_north 
    76    REAL, PUBLIC :: sign_north 
     85   REAL, PUBLIC    :: sign_north 
    7786   LOGICAL, PUBLIC :: l_ini_child = .FALSE. 
    78 # if defined key_vertical 
    79    LOGICAL, PUBLIC :: l_vremap    = .TRUE. 
    80 # else 
    8187   LOGICAL, PUBLIC :: l_vremap    = .FALSE. 
    82 # endif 
    8388!$AGRIF_END_DO_NOT_TREAT 
    8489    
    … …  
    100105      ALLOCATE( fspu(jpi,jpj), fspv(jpi,jpj),          & 
    101106         &      fspt(jpi,jpj), fspf(jpi,jpj),               & 
     107         &      fspu_2d(jpi,jpj), fspv_2d(jpi,jpj),         & 
     108         &      fspt_2d(jpi,jpj), fspf_2d(jpi,jpj),         & 
    102109         &      tabspongedone_tsn(jpi,jpj),                 & 
    103110         &      utint_stage(jpi,jpj), vtint_stage(jpi,jpj), & 
    … …  
    116123      ! 
    117124   END FUNCTION agrif_oce_alloc 
    118  
    119125#endif 
    120126   !!====================================================================== 

  • NEMO/trunk/src/NST/agrif_oce_interp.F90

    r14053 r14086  
    4545   PUBLIC   interpunb, interpvnb , interpub2b, interpvb2b 
    4646   PUBLIC   interpe3t, interpglamt, interpgphit 
    47    PUBLIC   interpht0, interpmbkt 
    48    PUBLIC   agrif_initts, agrif_initssh 
     47   PUBLIC   interpht0, interpmbkt, interpe3t0_vremap 
     48   PUBLIC   agrif_istate_oce, agrif_istate_ssh   ! called by icestate.F90 and domvvl.F90 
     49   PUBLIC   agrif_check_bat 
    4950 
    5051   INTEGER ::   bdy_tinterp = 0 
    … …  
    5859CONTAINS 
    5960 
    60    SUBROUTINE Agrif_tra 
    61       !!---------------------------------------------------------------------- 
    62       !!                  ***  ROUTINE Agrif_tra  *** 
    63       !!---------------------------------------------------------------------- 
    64       ! 
    65       IF( Agrif_Root() )   RETURN 
    66       ! 
     61   SUBROUTINE Agrif_istate_oce( Kbb, Kmm, Kaa ) 
     62      !!---------------------------------------------------------------------- 
     63      !!                 *** ROUTINE agrif_istate_oce *** 
     64      !! 
     65      !!                 set initial t, s, u, v, ssh from parent 
     66      !!---------------------------------------------------------------------- 
     67      ! 
     68      IMPLICIT NONE 
     69      ! 
     70      INTEGER, INTENT(in)  :: Kbb, Kmm, Kaa 
     71      INTEGER :: jn 
     72      !!---------------------------------------------------------------------- 
     73      IF(lwp) WRITE(numout,*) ' ' 
     74      IF(lwp) WRITE(numout,*) 'Agrif_istate_oce : interp child initial state from parent' 
     75      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~' 
     76      IF(lwp) WRITE(numout,*) ' ' 
     77 
     78      IF ( ln_rstart ) &  
     79         & CALL ctl_stop('AGRIF hot start should be desactivated in restarting mode') 
     80 
     81      IF ( .NOT.Agrif_Parent(l_1st_euler) ) &  
     82         & CALL ctl_stop('AGRIF hot start requires to force Euler first step on parent') 
     83 
     84      l_ini_child           = .TRUE. 
     85      Agrif_SpecialValue    = 0.0_wp 
     86      Agrif_UseSpecialValue = .TRUE. 
     87 
     88      ts(:,:,:,:,:) = 0.0_wp 
     89      uu(:,:,:,:)   = 0.0_wp 
     90      vv(:,:,:,:)   = 0.0_wp  
     91      ssh(:,:,:)    = 0._wp 
     92        
     93      Krhs_a = Kbb   ;   Kmm_a = Kbb 
     94 
     95      CALL Agrif_Init_Variable(tsini_id, procname=interptsn) 
     96      CALL Agrif_Init_Variable(sshini_id, procname=interpsshn) 
     97 
     98      Agrif_UseSpecialValue = ln_spc_dyn 
     99      use_sign_north = .TRUE. 
     100      sign_north = -1._wp 
     101      CALL Agrif_Init_Variable(uini_id , procname=interpun ) 
     102      CALL Agrif_Init_Variable(vini_id , procname=interpvn ) 
     103      use_sign_north = .FALSE. 
     104 
     105      Agrif_UseSpecialValue = .FALSE. 
     106      l_ini_child           = .FALSE. 
     107 
     108      Krhs_a = Kaa   ;   Kmm_a = Kmm 
     109 
     110      ssh(:,:,Kbb) = ssh(:,:,Kbb) * tmask(:,:,1) 
     111 
     112      DO jn = 1, jpts 
     113         ts(:,:,:,jn,Kbb) = ts(:,:,:,jn,Kbb) * tmask(:,:,:) 
     114      END DO 
     115      uu(:,:,:,Kbb) = uu(:,:,:,Kbb) * umask(:,:,:)      
     116      vv(:,:,:,Kbb) = vv(:,:,:,Kbb) * vmask(:,:,:)  
     117 
     118      CALL lbc_lnk_multi( 'agrif_istate_oce', uu(:,:,:  ,Kbb), 'U', -1.0_wp , vv(:,:,:,Kbb), 'V', -1.0_wp ) 
     119      CALL lbc_lnk( 'agrif_istate_oce', ts(:,:,:,:,Kbb), 'T',  1.0_wp ) 
     120      CALL lbc_lnk( 'agrif_istate_oce', ssh(:,:,Kbb), 'T', 1.0_wp ) 
     121 
     122   END SUBROUTINE Agrif_istate_oce 
     123 
     124 
     125   SUBROUTINE Agrif_istate_ssh( Kbb, Kmm ) 
     126      !!---------------------------------------------------------------------- 
     127      !!                 *** ROUTINE agrif_istate_ssh *** 
     128      !! 
     129      !!                    set initial ssh from parent 
     130      !!---------------------------------------------------------------------- 
     131      ! 
     132      IMPLICIT NONE 
     133      ! 
     134      INTEGER, INTENT(in)  :: Kbb, Kmm  
     135      !!---------------------------------------------------------------------- 
     136      IF(lwp) WRITE(numout,*) ' ' 
     137      IF(lwp) WRITE(numout,*) 'Agrif_istate_ssh : interp child ssh from parent' 
     138      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~' 
     139      IF(lwp) WRITE(numout,*) ' ' 
     140 
     141      IF ( ln_rstart ) &  
     142         & CALL ctl_stop('AGRIF hot start should be desactivated in restarting mode') 
     143 
     144      IF ( .NOT.Agrif_Parent(l_1st_euler) ) &  
     145         & CALL ctl_stop('AGRIF hot start requires to force Euler first step on parent') 
     146 
     147      Kmm_a = Kmm 
     148      ssh(:,:,Kmm) = 0._wp 
     149 
    67150      Agrif_SpecialValue    = 0._wp 
    68151      Agrif_UseSpecialValue = .TRUE. 
    69       ! 
    70       CALL Agrif_Bc_variable( tsn_id, procname=interptsn ) 
     152      l_ini_child           = .TRUE. 
     153      ! 
     154      CALL Agrif_Init_Variable(sshini_id, procname=interpsshn) 
    71155      ! 
    72156      Agrif_UseSpecialValue = .FALSE. 
     157      l_ini_child           = .FALSE. 
     158      CALL lbc_lnk( 'dom_vvl_rst', ssh(:,:,Kmm), 'T', 1._wp ) 
     159 
     160   END SUBROUTINE Agrif_istate_ssh 
     161 
     162 
     163   SUBROUTINE Agrif_tra 
     164      !!---------------------------------------------------------------------- 
     165      !!                  ***  ROUTINE Agrif_tra  *** 
     166      !!---------------------------------------------------------------------- 
     167      ! 
     168      IF( Agrif_Root() )   RETURN 
     169      ! 
     170      Agrif_SpecialValue    = 0._wp 
     171      Agrif_UseSpecialValue = .TRUE. 
     172      l_vremap           = ln_vert_remap 
     173      ! 
     174      CALL Agrif_Bc_variable( ts_interp_id, procname=interptsn ) 
     175      ! 
     176      Agrif_UseSpecialValue = .FALSE. 
     177      l_vremap              = .FALSE. 
    73178      ! 
    74179   END SUBROUTINE Agrif_tra 
    … …  
    90195      Agrif_SpecialValue    = 0.0_wp 
    91196      Agrif_UseSpecialValue = ln_spc_dyn 
     197      l_vremap              = ln_vert_remap 
    92198      ! 
    93199      use_sign_north = .TRUE. 
    … …  
    95201      CALL Agrif_Bc_variable( un_interp_id, procname=interpun ) 
    96202      CALL Agrif_Bc_variable( vn_interp_id, procname=interpvn ) 
     203 
     204      IF( .NOT.ln_dynspg_ts ) THEN ! Get transports 
     205         ubdy(:,:) = 0._wp   ;   vbdy(:,:) = 0._wp  
     206         CALL Agrif_Bc_variable( unb_interp_id, procname=interpunb ) 
     207         CALL Agrif_Bc_variable( vnb_interp_id, procname=interpvnb ) 
     208      ENDIF 
     209 
    97210      use_sign_north = .FALSE. 
    98211      ! 
    99212      Agrif_UseSpecialValue = .FALSE. 
     213      l_vremap              = .FALSE. 
     214      ! 
     215      ! Ensure below that vertically integrated transports match 
     216      ! either transports out of time splitting procedure (ln_dynspg_ts=.TRUE.) 
     217      ! or parent grid transports (ln_dynspg_ts=.FALSE.) 
    100218      ! 
    101219      ! --- West --- ! 
    102220      IF( lk_west ) THEN 
    103221         ibdy1 = nn_hls + 2                  ! halo + land + 1 
    104          ibdy2 = nn_hls + 1 + nbghostcells   ! halo + land + nbghostcells 
     222         ibdy2 = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhox()   ! halo + land + nbghostcells 
    105223         ! 
    106224         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport 
    107225            DO ji = mi0(ibdy1), mi1(ibdy2) 
    108                uu_b(ji,:,Krhs_a) = 0._wp 
    109                DO jk = 1, jpkm1 
    110                   DO jj = 1, jpj 
    111                      uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 
    112                   END DO 
    113                END DO 
    114226               DO jj = 1, jpj 
    115                   uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) * r1_hu(ji,jj,Krhs_a) 
     227                  uu_b(ji,jj,Krhs_a) = ubdy(ji,jj) * r1_hu(ji,jj,Krhs_a) 
     228                  vv_b(ji,jj,Krhs_a) = vbdy(ji,jj) * r1_hv(ji,jj,Krhs_a) 
    116229               END DO 
    117230            END DO 
    … …  
    119232         ! 
    120233         DO ji = mi0(ibdy1), mi1(ibdy2) 
    121             zub(ji,:) = 0._wp    ! Correct transport 
     234            zub(ji,:) = 0._wp   
    122235            DO jk = 1, jpkm1 
    123236               DO jj = 1, jpj 
    … …  
    135248         END DO 
    136249         !    
    137          IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate 
    138             DO ji = mi0(ibdy1), mi1(ibdy2) 
    139                zvb(ji,:) = 0._wp 
    140                DO jk = 1, jpkm1 
    141                   DO jj = 1, jpj 
    142                      zvb(ji,jj) = zvb(ji,jj) + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 
    143                   END DO 
    144                END DO 
     250         DO ji = mi0(ibdy1), mi1(ibdy2) 
     251            zvb(ji,:) = 0._wp 
     252            DO jk = 1, jpkm1 
    145253               DO jj = 1, jpj 
    146                   zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a) 
    147                END DO 
    148                DO jk = 1, jpkm1 
    149                   DO jj = 1, jpj 
    150                      vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) + vv_b(ji,jj,Krhs_a) - zvb(ji,jj) )*vmask(ji,jj,jk) 
    151                   END DO 
    152                END DO 
    153             END DO 
    154          ENDIF 
     254                  zvb(ji,jj) = zvb(ji,jj) + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 
     255               END DO 
     256            END DO 
     257            DO jj = 1, jpj 
     258               zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a) 
     259            END DO 
     260            DO jk = 1, jpkm1 
     261               DO jj = 1, jpj 
     262                  vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) + vv_b(ji,jj,Krhs_a) - zvb(ji,jj) )*vmask(ji,jj,jk) 
     263               END DO 
     264            END DO 
     265         END DO 
    155266         ! 
    156267      ENDIF 
    … …  
    158269      ! --- East --- ! 
    159270      IF( lk_east) THEN 
    160          ibdy1 = jpiglo - ( nn_hls + nbghostcells + 1)   ! halo + land + nbghostcells 
    161          ibdy2 = jpiglo - ( nn_hls + 2 )                 ! halo + land + 1 
    162          ! 
    163          IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport 
     271         ibdy1 = jpiglo - ( nn_hls + nbghostcells + 1) - nn_shift_bar*Agrif_Rhox()     
     272         ibdy2 = jpiglo - ( nn_hls + 2 )                  
     273         ! 
     274         IF( .NOT.ln_dynspg_ts ) THEN  
    164275            DO ji = mi0(ibdy1), mi1(ibdy2) 
    165276               uu_b(ji,:,Krhs_a) = 0._wp 
    … …  
    170281               END DO 
    171282               DO jj = 1, jpj 
    172                   uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) * r1_hu(ji,jj,Krhs_a) 
     283                  uu_b(ji,jj,Krhs_a) = ubdy(ji,jj) * r1_hu(ji,jj,Krhs_a) 
    173284               END DO 
    174285            END DO 
    … …  
    176287         ! 
    177288         DO ji = mi0(ibdy1), mi1(ibdy2) 
    178             zub(ji,:) = 0._wp    ! Correct transport 
     289            zub(ji,:) = 0._wp    
    179290            DO jk = 1, jpkm1 
    180291               DO jj = 1, jpj 
    … …  
    192303         END DO 
    193304         ! 
    194          IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate 
    195             ibdy1 = jpiglo - ( nn_hls + nbghostcells )   ! halo + land + nbghostcells - 1 
    196             ibdy2 = jpiglo - ( nn_hls + 1 )              ! halo + land + 1            - 1 
     305         ibdy1 = jpiglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhox()  
     306         ibdy2 = jpiglo - ( nn_hls + 1 )              ! 
     307         IF( .NOT.ln_dynspg_ts ) THEN  
    197308            DO ji = mi0(ibdy1), mi1(ibdy2) 
    198                zvb(ji,:) = 0._wp 
     309               vv_b(ji,:,Krhs_a) = 0._wp 
    199310               DO jk = 1, jpkm1 
    200311                  DO jj = 1, jpj 
     312                     vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 
     313                  END DO 
     314               END DO 
     315               DO jj = 1, jpj 
     316                  vv_b(ji,jj,Krhs_a) = vbdy(ji,jj) * r1_hv(ji,jj,Krhs_a) 
     317               END DO 
     318            END DO 
     319         ENDIF 
     320 
     321         DO ji = mi0(ibdy1), mi1(ibdy2) 
     322            zvb(ji,:) = 0._wp 
     323            DO jk = 1, jpkm1 
     324               DO jj = 1, jpj 
    201325                     zvb(ji,jj) = zvb(ji,jj) + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 
    202                   END DO 
    203                END DO 
     326               END DO 
     327            END DO 
     328            DO jj = 1, jpj 
     329               zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a) 
     330            END DO 
     331            DO jk = 1, jpkm1 
    204332               DO jj = 1, jpj 
    205                   zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a) 
    206                END DO 
    207                DO jk = 1, jpkm1 
    208                   DO jj = 1, jpj 
    209333                     vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) + vv_b(ji,jj,Krhs_a) - zvb(ji,jj) ) * vmask(ji,jj,jk) 
    210                   END DO 
    211                END DO 
    212             END DO 
    213          ENDIF 
     334               END DO 
     335            END DO 
     336         END DO 
    214337         ! 
    215338      ENDIF 
    … …  
    217340      ! --- South --- ! 
    218341      IF( lk_south ) THEN 
    219          jbdy1 = nn_hls + 2                  ! halo + land + 1 
    220          jbdy2 = nn_hls + 1 + nbghostcells   ! halo + land + nbghostcells 
    221          ! 
    222          IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport 
     342         jbdy1 = nn_hls + 2                  
     343         jbdy2 = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhoy()    
     344         ! 
     345         IF( .NOT.ln_dynspg_ts ) THEN 
    223346            DO jj = mj0(jbdy1), mj1(jbdy2) 
    224347               vv_b(:,jj,Krhs_a) = 0._wp 
    … …  
    229352               END DO 
    230353               DO ji=1,jpi 
    231                   vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) * r1_hv(ji,jj,Krhs_a)      
     354                  vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) * r1_hv(ji,jj,Krhs_a)  
     355                  uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) * r1_hu(ji,jj,Krhs_a) 
    232356               END DO 
    233357            END DO 
    … …  
    235359         ! 
    236360         DO jj = mj0(jbdy1), mj1(jbdy2) 
    237             zvb(:,jj) = 0._wp    ! Correct transport 
     361            zvb(:,jj) = 0._wp 
    238362            DO jk=1,jpkm1 
    239363               DO ji=1,jpi 
    … …  
    251375         END DO 
    252376         ! 
    253          IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate 
    254             DO jj = mj0(jbdy1), mj1(jbdy2) 
    255                zub(:,jj) = 0._wp 
    256                DO jk = 1, jpkm1 
    257                   DO ji = 1, jpi 
    258                      zub(ji,jj) = zub(ji,jj) + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 
    259                   END DO 
    260                END DO 
     377         DO jj = mj0(jbdy1), mj1(jbdy2) 
     378            zub(:,jj) = 0._wp 
     379            DO jk = 1, jpkm1 
    261380               DO ji = 1, jpi 
    262                   zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 
    263                END DO 
    264                DO jk = 1, jpkm1 
    265                   DO ji = 1, jpi 
    266                      uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) + uu_b(ji,jj,Krhs_a) - zub(ji,jj) ) * umask(ji,jj,jk) 
    267                   END DO 
    268                END DO 
    269             END DO 
    270          ENDIF 
     381                  zub(ji,jj) = zub(ji,jj) + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 
     382               END DO 
     383            END DO 
     384            DO ji = 1, jpi 
     385               zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 
     386            END DO 
     387            DO jk = 1, jpkm1 
     388               DO ji = 1, jpi 
     389                  uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) + uu_b(ji,jj,Krhs_a) - zub(ji,jj) ) * umask(ji,jj,jk) 
     390               END DO 
     391            END DO 
     392         END DO 
    271393         ! 
    272394      ENDIF 
    … …  
    274396      ! --- North --- ! 
    275397      IF( lk_north ) THEN 
    276          jbdy1 = jpjglo - ( nn_hls + nbghostcells + 1)   ! halo + land + nbghostcells 
    277          jbdy2 = jpjglo - ( nn_hls + 2 )                 ! halo + land + 1 
    278          ! 
    279          IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport 
     398         jbdy1 = jpjglo - ( nn_hls + nbghostcells + 1) - nn_shift_bar*Agrif_Rhoy()  
     399         jbdy2 = jpjglo - ( nn_hls + 2 ) 
     400         ! 
     401         IF( .NOT.ln_dynspg_ts ) THEN 
    280402            DO jj = mj0(jbdy1), mj1(jbdy2) 
    281403               vv_b(:,jj,Krhs_a) = 0._wp 
    … …  
    286408               END DO 
    287409               DO ji=1,jpi 
    288                   vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) * r1_hv(ji,jj,Krhs_a) 
     410                  vv_b(ji,jj,Krhs_a) = vbdy(ji,jj) * r1_hv(ji,jj,Krhs_a) 
    289411               END DO 
    290412            END DO 
    … …  
    292414         ! 
    293415         DO jj = mj0(jbdy1), mj1(jbdy2) 
    294             zvb(:,jj) = 0._wp    ! Correct transport 
     416            zvb(:,jj) = 0._wp  
    295417            DO jk=1,jpkm1 
    296418               DO ji=1,jpi 
    … …  
    308430         END DO 
    309431         ! 
    310          IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate 
    311             jbdy1 = jpjglo - ( nn_hls + nbghostcells )   ! halo + land + nbghostcells - 1 
    312             jbdy2 = jpjglo - ( nn_hls + 1 )              ! halo + land + 1            - 1 
     432         jbdy1 = jpjglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhoy()   
     433         jbdy2 = jpjglo - ( nn_hls + 1 ) 
     434         IF( .NOT.ln_dynspg_ts ) THEN 
    313435            DO jj = mj0(jbdy1), mj1(jbdy2) 
    314                zub(:,jj) = 0._wp 
     436               uu_b(:,jj,Krhs_a) = 0._wp 
    315437               DO jk = 1, jpkm1 
    316438                  DO ji = 1, jpi 
    317                      zub(ji,jj) = zub(ji,jj) + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 
    318                   END DO 
    319                END DO 
     439                     uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 
     440                  END DO 
     441               END DO 
     442               DO ji=1,jpi 
     443                  uu_b(ji,jj,Krhs_a) = ubdy(ji,jj) * r1_hu(ji,jj,Krhs_a) 
     444               END DO 
     445            END DO 
     446         ENDIF 
     447 
     448         DO jj = mj0(jbdy1), mj1(jbdy2) 
     449            zub(:,jj) = 0._wp 
     450            DO jk = 1, jpkm1 
    320451               DO ji = 1, jpi 
    321                   zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 
    322                END DO 
    323                DO jk = 1, jpkm1 
    324                   DO ji = 1, jpi 
     452                  zub(ji,jj) = zub(ji,jj) + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 
     453               END DO 
     454            END DO 
     455            DO ji = 1, jpi 
     456               zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 
     457            END DO 
     458            DO jk = 1, jpkm1 
     459               DO ji = 1, jpi 
    325460                     uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) + uu_b(ji,jj,Krhs_a) - zub(ji,jj) ) * umask(ji,jj,jk) 
    326                   END DO 
    327                END DO 
    328             END DO 
    329          ENDIF 
     461               END DO 
     462            END DO 
     463         END DO 
    330464         ! 
    331465      ENDIF 
    … …  
    349483      IF( lk_west ) THEN 
    350484         istart = nn_hls + 2                              ! halo + land + 1 
    351          iend   = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
     485         iend   = nn_hls + 1 + nbghostcells  + nn_shift_bar*Agrif_Rhox()              ! halo + land + nbghostcells 
    352486         DO ji = mi0(istart), mi1(iend) 
    353487            DO jj=1,jpj 
    … …  
    360494      !--- East ---! 
    361495      IF( lk_east ) THEN 
    362          istart = jpiglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    363          iend   = jpiglo - ( nn_hls + 1 )                 ! halo + land + 1            - 1 
     496         istart = jpiglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhox()  
     497         iend   = jpiglo - ( nn_hls + 1 )                 
    364498         DO ji = mi0(istart), mi1(iend) 
    365499 
    … …  
    368502            END DO 
    369503         END DO 
    370          istart = jpiglo - ( nn_hls + nbghostcells + 1)   ! halo + land + nbghostcells 
    371          iend   = jpiglo - ( nn_hls + 2 )                 ! halo + land + 1 
     504         istart = jpiglo - ( nn_hls + nbghostcells + 1) - nn_shift_bar*Agrif_Rhox()  
     505         iend   = jpiglo - ( nn_hls + 2 )                 
    372506         DO ji = mi0(istart), mi1(iend) 
    373507            DO jj=1,jpj 
    … …  
    379513      !--- South ---! 
    380514      IF( lk_south ) THEN 
    381          jstart = nn_hls + 2                              ! halo + land + 1 
    382          jend   = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
     515         jstart = nn_hls + 2                               
     516         jend   = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhoy()            
    383517         DO jj = mj0(jstart), mj1(jend) 
    384518 
    … …  
    392526      !--- North ---! 
    393527      IF( lk_north ) THEN 
    394          jstart = jpjglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    395          jend   = jpjglo - ( nn_hls + 1 )                 ! halo + land + 1            - 1 
     528         jstart = jpjglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhoy()      
     529         jend   = jpjglo - ( nn_hls + 1 )                 
    396530         DO jj = mj0(jstart), mj1(jend) 
    397531            DO ji=1,jpi 
    … …  
    399533            END DO 
    400534         END DO 
    401          jstart = jpjglo - ( nn_hls + nbghostcells + 1)   ! halo + land + nbghostcells 
    402          jend   = jpjglo - ( nn_hls + 2 )                 ! halo + land + 1 
     535         jstart = jpjglo - ( nn_hls + nbghostcells + 1) - nn_shift_bar*Agrif_Rhoy()  
     536         jend   = jpjglo - ( nn_hls + 2 )                 
    403537         DO jj = mj0(jstart), mj1(jend) 
    404538            DO ji=1,jpi 
    … …  
    426560      !--- West ---! 
    427561      IF( lk_west ) THEN 
    428          istart = nn_hls + 2                              ! halo + land + 1 
    429          iend   = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
     562         istart = nn_hls + 2                               
     563         iend   = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhox()  
    430564         DO ji = mi0(istart), mi1(iend) 
    431565            DO jj=1,jpj 
    … …  
    438572      !--- East ---! 
    439573      IF( lk_east ) THEN 
    440          istart = jpiglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    441          iend   = jpiglo - ( nn_hls + 1 )                 ! halo + land + 1            - 1 
     574         istart = jpiglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhox() 
     575         iend   = jpiglo - ( nn_hls + 1 )                  
    442576         DO ji = mi0(istart), mi1(iend) 
    443577            DO jj=1,jpj 
    … …  
    445579            END DO 
    446580         END DO 
    447          istart = jpiglo - ( nn_hls + nbghostcells + 1)   ! halo + land + nbghostcells 
    448          iend   = jpiglo - ( nn_hls + 2 )                 ! halo + land + 1 
     581         istart = jpiglo - ( nn_hls + nbghostcells + 1) - nn_shift_bar*Agrif_Rhox()  
     582         iend   = jpiglo - ( nn_hls + 2 )                  
    449583         DO ji = mi0(istart), mi1(iend) 
    450584            DO jj=1,jpj 
    … …  
    456590      !--- South ---! 
    457591      IF( lk_south ) THEN 
    458          jstart = nn_hls + 2                              ! halo + land + 1 
    459          jend   = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
     592         jstart = nn_hls + 2                               
     593         jend   = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhoy()  
    460594         DO jj = mj0(jstart), mj1(jend) 
    461595            DO ji=1,jpi 
    … …  
    468602      !--- North ---! 
    469603      IF( lk_north ) THEN 
    470          jstart = jpjglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    471          jend   = jpjglo - ( nn_hls + 1 )                 ! halo + land + 1            - 1 
     604         jstart = jpjglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhoy()  
     605         jend   = jpjglo - ( nn_hls + 1 )                 
    472606         DO jj = mj0(jstart), mj1(jend) 
    473607            DO ji=1,jpi 
    … …  
    475609            END DO 
    476610         END DO 
    477          jstart = jpjglo - ( nn_hls + nbghostcells + 1)   ! halo + land + nbghostcells 
    478          jend   = jpjglo - ( nn_hls + 2 )                 ! halo + land + 1 
     611         jstart = jpjglo - ( nn_hls + nbghostcells + 1) - nn_shift_bar*Agrif_Rhoy()  
     612         jend   = jpjglo - ( nn_hls + 2 )                
    479613         DO jj = mj0(jstart), mj1(jend) 
    480614            DO ji=1,jpi 
    … …  
    493627      INTEGER, INTENT(in) ::   kt 
    494628      !! 
    495       INTEGER :: ji, jj 
    496629      LOGICAL :: ll_int_cons 
    497630      !!----------------------------------------------------------------------   
    … …  
    517650      ! 
    518651      IF( ll_int_cons ) THEN  ! Conservative interpolation 
    519          ! order matters here !!!!!! 
    520          CALL Agrif_Bc_variable( ub2b_interp_id, calledweight=1._wp, procname=interpub2b ) ! Time integrated 
    521          CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b )  
    522          ! 
    523          bdy_tinterp = 1 
    524          CALL Agrif_Bc_variable( unb_id        , calledweight=1._wp, procname=interpunb  ) ! After 
    525          CALL Agrif_Bc_variable( vnb_id        , calledweight=1._wp, procname=interpvnb  )   
    526          ! 
    527          bdy_tinterp = 2 
    528          CALL Agrif_Bc_variable( unb_id        , calledweight=0._wp, procname=interpunb  ) ! Before 
    529          CALL Agrif_Bc_variable( vnb_id        , calledweight=0._wp, procname=interpvnb  )    
     652         IF ( lk_tint2d_notinterp ) THEN 
     653            CALL Agrif_Bc_variable( ub2b_interp_id, calledweight=1._wp, procname=interpub2b_const ) 
     654            CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b_const )  
     655            ! Divergence conserving correction terms: 
     656            CALL Agrif_Bc_variable(    ub2b_cor_id, calledweight=1._wp, procname=        ub2b_cor ) 
     657            CALL Agrif_Bc_variable(    vb2b_cor_id, calledweight=1._wp, procname=        vb2b_cor ) 
     658         ELSE 
     659            ! order matters here !!!!!! 
     660            CALL Agrif_Bc_variable( ub2b_interp_id, calledweight=1._wp, procname=interpub2b ) ! Time integrated 
     661            CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b )  
     662            ! 
     663            bdy_tinterp = 1 
     664            CALL Agrif_Bc_variable( unb_interp_id , calledweight=1._wp, procname=interpunb  ) ! After 
     665            CALL Agrif_Bc_variable( vnb_interp_id , calledweight=1._wp, procname=interpvnb  )   
     666            ! 
     667            bdy_tinterp = 2 
     668            CALL Agrif_Bc_variable( unb_interp_id , calledweight=0._wp, procname=interpunb  ) ! Before 
     669            CALL Agrif_Bc_variable( vnb_interp_id , calledweight=0._wp, procname=interpvnb  )    
     670         ENDIF 
    530671      ELSE ! Linear interpolation 
    531672         ! 
    532673         ubdy(:,:) = 0._wp   ;   vbdy(:,:) = 0._wp  
    533          CALL Agrif_Bc_variable( unb_id, procname=interpunb ) 
    534          CALL Agrif_Bc_variable( vnb_id, procname=interpvnb ) 
     674         CALL Agrif_Bc_variable( unb_interp_id, procname=interpunb ) 
     675         CALL Agrif_Bc_variable( vnb_interp_id, procname=interpvnb ) 
    535676      ENDIF 
    536677      Agrif_UseSpecialValue = .FALSE. 
    … …  
    561702      ! --- West --- ! 
    562703      IF(lk_west) THEN 
    563          istart = nn_hls + 2                              ! halo + land + 1 
    564          iend   = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
     704         istart = nn_hls + 2                                                          ! halo + land + 1 
     705         iend   = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhox()               ! halo + land + nbghostcells 
    565706         DO ji = mi0(istart), mi1(iend) 
    566707            DO jj = 1, jpj 
    … …  
    572713      ! --- East --- ! 
    573714      IF(lk_east) THEN 
    574          istart = jpiglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    575          iend   = jpiglo - ( nn_hls + 1 )                 ! halo + land + 1            - 1 
     715         istart = jpiglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhox()       ! halo + land + nbghostcells - 1 
     716         iend   = jpiglo - ( nn_hls + 1 )                                              ! halo + land + 1            - 1 
    576717         DO ji = mi0(istart), mi1(iend) 
    577718            DO jj = 1, jpj 
    … …  
    583724      ! --- South --- ! 
    584725      IF(lk_south) THEN 
    585          jstart = nn_hls + 2                              ! halo + land + 1 
    586          jend   = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
     726         jstart = nn_hls + 2                                                          ! halo + land + 1 
     727         jend   = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhoy()               ! halo + land + nbghostcells 
    587728         DO jj = mj0(jstart), mj1(jend) 
    588729            DO ji = 1, jpi 
    … …  
    594735      ! --- North --- ! 
    595736      IF(lk_north) THEN 
    596          jstart = jpjglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    597          jend   = jpjglo - ( nn_hls + 1 )                 ! halo + land + 1            - 1 
     737         jstart = jpjglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhoy()     ! halo + land + nbghostcells - 1 
     738         jend   = jpjglo - ( nn_hls + 1 )                                            ! halo + land + 1            - 1 
    598739         DO jj = mj0(jstart), mj1(jend) 
    599740            DO ji = 1, jpi 
    … …  
    620761      ! --- West --- ! 
    621762      IF(lk_west) THEN 
    622          istart = nn_hls + 2                              ! halo + land + 1 
    623          iend   = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
     763         istart = nn_hls + 2                                                        ! halo + land + 1 
     764         iend   = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhox()             ! halo + land + nbghostcells 
    624765         DO ji = mi0(istart), mi1(iend) 
    625766            DO jj = 1, jpj 
    … …  
    631772      ! --- East --- ! 
    632773      IF(lk_east) THEN 
    633          istart = jpiglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    634          iend   = jpiglo - ( nn_hls + 1 )                 ! halo + land + 1            - 1 
     774         istart = jpiglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhox()    ! halo + land + nbghostcells - 1 
     775         iend   = jpiglo - ( nn_hls + 1 )                                           ! halo + land + 1            - 1 
    635776         DO ji = mi0(istart), mi1(iend) 
    636777            DO jj = 1, jpj 
    … …  
    642783      ! --- South --- ! 
    643784      IF(lk_south) THEN 
    644          jstart = nn_hls + 2                              ! halo + land + 1 
    645          jend   = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
     785         jstart = nn_hls + 2                                                        ! halo + land + 1 
     786         jend   = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhoy()             ! halo + land + nbghostcells 
    646787         DO jj = mj0(jstart), mj1(jend) 
    647788            DO ji = 1, jpi 
    … …  
    653794      ! --- North --- ! 
    654795      IF(lk_north) THEN 
    655          jstart = jpjglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    656          jend   = jpjglo - ( nn_hls + 1 )                 ! halo + land + 1            - 1 
     796         jstart = jpjglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhoy()    ! halo + land + nbghostcells - 1 
     797         jend   = jpjglo - ( nn_hls + 1 )                                           ! halo + land + 1            - 1 
    657798         DO jj = mj0(jstart), mj1(jend) 
    658799            DO ji = 1, jpi 
    … …  
    679820      Agrif_SpecialValue    = 0.e0 
    680821      Agrif_UseSpecialValue = .TRUE. 
     822      l_vremap              = ln_vert_remap 
    681823      ! 
    682824      CALL Agrif_Bc_variable( avm_id, calledweight=zalpha, procname=interpavm )        
    683825      ! 
    684826      Agrif_UseSpecialValue = .FALSE. 
     827      l_vremap              = .FALSE. 
    685828      ! 
    686829   END SUBROUTINE Agrif_avm 
    … …  
    688831 
    689832   SUBROUTINE interptsn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    690       !!---------------------------------------------------------------------- 
    691       !!                  *** ROUTINE interptsn *** 
    692833      !!---------------------------------------------------------------------- 
    693834      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   ptab 
    … …  
    699840      INTEGER  :: item 
    700841      ! vertical interpolation: 
    701       REAL(wp) :: zhtot 
    702       REAL(wp), DIMENSION(k1:k2,1:jpts) :: tabin 
    703       REAL(wp), DIMENSION(k1:k2) :: h_in, z_in 
     842      REAL(wp) :: zhtot, zwgt 
     843      REAL(wp), DIMENSION(k1:k2,1:jpts) :: tabin, tabin_i 
     844      REAL(wp), DIMENSION(k1:k2) :: z_in, h_in_i, z_in_i 
    704845      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out 
    705846      !!---------------------------------------------------------------------- 
    … …  
    720861         END DO 
    721862 
    722          IF( l_vremap .OR. l_ini_child) THEN 
    723             ! Interpolate thicknesses 
     863         IF( l_vremap .OR. l_ini_child .OR. ln_zps ) THEN 
     864 
     865            ! Fill cell depths (i.e. gdept) to be interpolated 
    724866            ! Warning: these are masked, hence extrapolated prior interpolation. 
    725             DO jk=k1,k2 
    726                DO jj=j1,j2 
    727                   DO ji=i1,i2 
    728                       ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 
    729  
    730                   END DO 
    731                END DO 
    732             END DO 
    733  
    734             ! Extrapolate thicknesses in partial bottom cells: 
    735             ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
    736             IF (ln_zps) THEN 
    737                DO jj=j1,j2 
    738                   DO ji=i1,i2 
    739                       jk = mbkt(ji,jj) 
    740                       ptab(ji,jj,jk,jpts+1) = 0._wp 
    741                   END DO 
    742                END DO            
    743             END IF 
    744          
     867            DO jj=j1,j2 
     868               DO ji=i1,i2 
     869                  ptab(ji,jj,k1,jpts+1) = 0.5_wp * tmask(ji,jj,k1) * e3t(ji,jj,k1,Kmm_a) 
     870                  DO jk=k1+1,k2 
     871                     ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * & 
     872                        & ( ptab(ji,jj,jk-1,jpts+1) + 0.5_wp * (e3t(ji,jj,jk-1,Kmm_a)+e3t(ji,jj,jk,Kmm_a)) ) 
     873                  END DO 
     874               END DO 
     875            END DO 
     876          
    745877            ! Save ssh at last level: 
    746878            IF (.NOT.ln_linssh) THEN 
    747879               ptab(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1)  
    748             ELSE 
    749                ptab(i1:i2,j1:j2,k2,jpts+1) = 0._wp 
    750880            END IF       
    751881         ENDIF 
    … …  
    758888         IF( l_vremap .OR. l_ini_child ) THEN 
    759889            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,n2) = 0._wp  
    760                 
    761890            DO jj=j1,j2 
    762891               DO ji=i1,i2 
    763                   ts(ji,jj,:,:,Krhs_a) = 0.                   
    764                !   IF( l_ini_child) ts(ji,jj,:,:,Krhs_a) = ptab(ji,jj,:,1:jpts) 
     892                  ts(ji,jj,:,:,Krhs_a) = 0.   
     893                  ! 
     894                  ! Build vertical grids: 
    765895                  N_in = mbkt_parent(ji,jj) 
    766                   zhtot = 0._wp 
    767                   DO jk=1,N_in !k2 = jpk of parent grid 
    768                      IF (jk==N_in) THEN 
    769                         h_in(jk) = ht0_parent(ji,jj) + ptab(ji,jj,k2,n2) - zhtot 
    770                      ELSE 
    771                         h_in(jk) = ptab(ji,jj,jk,n2) 
    772                      ENDIF 
    773                      zhtot = zhtot + h_in(jk) 
    774                      tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1) 
    775                   END DO 
    776                   z_in(1) = 0.5_wp * h_in(1) - zhtot + ht0_parent(ji,jj) 
    777                   DO jk=2,N_in 
    778                      z_in(jk) = z_in(jk-1) + 0.5_wp * h_in(jk) 
    779                   END DO 
    780  
    781                   N_out = 0 
    782                   DO jk=1,jpk ! jpk of child grid 
    783                      IF (tmask(ji,jj,jk) == 0._wp) EXIT  
    784                      N_out = N_out + 1 
     896                  ! Input grid (account for partial cells if any): 
     897                  DO jk=1,N_in 
     898                     z_in(jk) = ptab(ji,jj,jk,n2) - ptab(ji,jj,k2,n2) 
     899                     tabin(jk,1:jpts) = ptab(ji,jj,jk,1:jpts) 
     900                  END DO 
     901                   
     902                  ! Intermediate grid: 
     903                  IF ( l_vremap ) THEN 
     904                     DO jk = 1, N_in 
     905                        h_in_i(jk) = e3t0_parent(ji,jj,jk) * &  
     906                          &       (1._wp + ptab(ji,jj,k2,n2)/(ht0_parent(ji,jj)*ssmask(ji,jj) + 1._wp - ssmask(ji,jj))) 
     907                     END DO 
     908                     z_in_i(1) = 0.5_wp * h_in_i(1) 
     909                     DO jk=2,N_in 
     910                        z_in_i(jk) = z_in_i(jk-1) + 0.5_wp * ( h_in_i(jk) + h_in_i(jk-1) ) 
     911                     END DO 
     912                     z_in_i(1:N_in) = z_in_i(1:N_in)  - ptab(ji,jj,k2,n2) 
     913                  ENDIF                               
     914 
     915                  ! Output (Child) grid: 
     916                  N_out = mbkt(ji,jj) 
     917                  DO jk=1,N_out 
    785918                     h_out(jk) = e3t(ji,jj,jk,Krhs_a) 
    786919                  END DO 
    787  
    788                   z_out(1) = 0.5_wp * h_out(1) - SUM(h_out(1:N_out)) + ht_0(ji,jj) 
     920                  z_out(1) = 0.5_wp * h_out(1) 
    789921                  DO jk=2,N_out 
    790                      z_out(jk) = z_out(jk-1) + 0.5_wp * h_out(jk) 
    791                   END DO 
     922                     z_out(jk) = z_out(jk-1) + 0.5_wp * ( h_out(jk)+h_out(jk-1) ) 
     923                  END DO 
     924                  IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out)  - ssh(ji,jj,Krhs_a) 
    792925 
    793926                  IF (N_in*N_out > 0) THEN 
    794927                     IF( l_ini_child ) THEN 
    795                         CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),          & 
     928                        CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),              & 
    796929                                      &   z_out(1:N_out),N_in,N_out,jpts)   
    797930                     ELSE  
    798                         CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),   & 
     931                        CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),tabin_i(1:N_in,1:jpts),                       & 
     932                                     &   z_in_i(1:N_in),N_in,N_in,jpts) 
     933                        CALL reconstructandremap(tabin_i(1:N_in,1:jpts),h_in_i(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),   & 
    799934                                      &   h_out(1:N_out),N_in,N_out,jpts)   
    800935                     ENDIF 
    … …  
    806941         ELSE 
    807942          
    808             DO jn=1, jpts 
    809                 ts(i1:i2,j1:j2,1:jpk,jn,Krhs_a)=ptab(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk)  
     943            IF ( Agrif_Parent(ln_zps) ) THEN ! Account for partial cells  
     944                                             ! linear vertical interpolation 
     945               DO jj=j1,j2 
     946                  DO ji=i1,i2 
     947                     ! 
     948                     N_in  = mbkt(ji,jj) 
     949                     N_out = mbkt(ji,jj) 
     950                     z_in(1) = ptab(ji,jj,1,n2) 
     951                     tabin(1,1:jpts) = ptab(ji,jj,1,1:jpts) 
     952                     DO jk=2, N_in 
     953                        z_in(jk) = ptab(ji,jj,jk,n2) 
     954                        tabin(jk,1:jpts) = ptab(ji,jj,jk,1:jpts) 
     955                     END DO 
     956                     IF (.NOT.ln_linssh) z_in(1:N_in) = z_in(1:N_in) - ptab(ji,jj,k2,n2) 
     957                     z_out(1) = 0.5_wp * e3t(ji,jj,1,Krhs_a) 
     958                     DO jk=2, N_out 
     959                        z_out(jk) = z_out(jk-1) + 0.5_wp * (e3t(ji,jj,jk-1,Krhs_a) + e3t(ji,jj,jk,Krhs_a)) 
     960                     END DO 
     961                     IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out) - ssh(ji,jj,Krhs_a) 
     962                     CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),ptab(ji,jj,1:N_out,1:jpts), & 
     963                                   &   z_out(1:N_out),N_in,N_out,jpts)   
     964                  END DO 
     965               END DO 
     966            ENDIF 
     967 
     968            DO jn =1, jpts 
     969               ts(i1:i2,j1:j2,1:jpk,jn,Krhs_a) = ptab(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 
    810970            END DO 
    811971         ENDIF 
    … …  
    829989         ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Kmm_a) 
    830990      ELSE 
    831          hbdy(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1) 
     991         IF( l_ini_child ) THEN 
     992            ssh(i1:i2,j1:j2,Kmm_a) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1) 
     993         ELSE 
     994            hbdy(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1) 
     995         ENDIF 
    832996      ENDIF 
    833997      ! 
    … …  
    8701034         END DO 
    8711035 
    872         IF( l_vremap .OR. l_ini_child) THEN 
     1036        IF( l_vremap .OR. l_ini_child ) THEN 
    8731037         ! Extrapolate thicknesses in partial bottom cells: 
    8741038         ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
    … …  
    9091073                  zhtot = 0._wp 
    9101074                  DO jk=1,N_in 
    911                      IF (jk==N_in) THEN 
    912                         h_in(jk) = hu0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot 
     1075                     !IF (jk==N_in) THEN 
     1076                     !   h_in(jk) = hu0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot 
     1077                     !ELSE 
     1078                     !   h_in(jk) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy)  
     1079                     !ENDIF 
     1080                     IF ( l_vremap ) THEN 
     1081                        h_in(jk) = e3u0_parent(ji,jj,jk)  
    9131082                     ELSE 
    914                         h_in(jk) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy)  
     1083                        IF (jk==N_in) THEN 
     1084                           h_in(jk) = hu0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot 
     1085                        ELSE 
     1086                           h_in(jk) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy)  
     1087                        ENDIF 
    9151088                     ENDIF 
    9161089                     zhtot = zhtot + h_in(jk) 
    … …  
    9231096                 z_in(1) = 0.5_wp * h_in(1) - zhtot + hu0_parent(ji,jj)  
    9241097                 DO jk=2,N_in 
    925                     z_in(jk) = z_in(jk-1) + 0.5_wp * h_in(jk) 
     1098                    z_in(jk) = z_in(jk-1) + 0.5_wp * (h_in(jk)+h_in(jk-1)) 
    9261099                 END DO 
    9271100                      
    … …  
    9351108                 z_out(1) = 0.5_wp * h_out(1) - SUM(h_out(1:N_out)) + hu_0(ji,jj) 
    9361109                 DO jk=2,N_out 
    937                     z_out(jk) = z_out(jk-1) + 0.5_wp * h_out(jk)  
     1110                    z_out(jk) = z_out(jk-1) + 0.5_wp * (h_out(jk-1) + h_out(jk))  
    9381111                 END DO   
    9391112 
    … …  
    10311204                  zhtot = 0._wp 
    10321205                  DO jk=1,N_in 
    1033                      IF (jk==N_in) THEN 
    1034                         h_in(jk) = hv0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot 
     1206                     !IF (jk==N_in) THEN 
     1207                     !   h_in(jk) = hv0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot 
     1208                     !ELSE 
     1209                     !   h_in(jk) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox)  
     1210                     !ENDIF 
     1211                     IF (l_vremap) THEN 
     1212                        h_in(jk) = e3v0_parent(ji,jj,jk) 
    10351213                     ELSE 
    1036                         h_in(jk) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox)  
     1214                        IF (jk==N_in) THEN 
     1215                           h_in(jk) = hv0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot 
     1216                        ELSE 
     1217                           h_in(jk) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox)  
     1218                        ENDIF 
    10371219                     ENDIF 
    10381220                     zhtot = zhtot + h_in(jk) 
    … …  
    10461228                  z_in(1) = 0.5_wp * h_in(1) - zhtot + hv0_parent(ji,jj) 
    10471229                  DO jk=2,N_in 
    1048                      z_in(jk) = z_in(jk-1) + 0.5_wp * h_in(jk) 
     1230                     z_in(jk) = z_in(jk-1) + 0.5_wp * (h_in(jk-1)+h_in(jk)) 
    10491231                  END DO 
    10501232 
    … …  
    10581240                  z_out(1) = 0.5_wp * h_out(1) - SUM(h_out(1:N_out)) + hv_0(ji,jj) 
    10591241                  DO jk=2,N_out 
    1060                      z_out(jk) = z_out(jk-1) + 0.5_wp * h_out(jk) 
     1242                     z_out(jk) = z_out(jk-1) + 0.5_wp * (h_out(jk-1)+h_out(jk)) 
    10611243                  END DO 
    10621244  
    … …  
    12151397   END SUBROUTINE interpub2b 
    12161398    
     1399   SUBROUTINE interpub2b_const( ptab, i1, i2, j1, j2, before ) 
     1400      !!---------------------------------------------------------------------- 
     1401      !!                  ***  ROUTINE interpub2b_const  *** 
     1402      !!----------------------------------------------------------------------   
     1403      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     1404      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     1405      LOGICAL                         , INTENT(in   ) ::   before 
     1406      ! 
     1407      REAL(wp) :: zrhoy 
     1408      !!----------------------------------------------------------------------   
     1409      IF( before ) THEN 
     1410!         IF ( ln_bt_fw ) THEN 
     1411            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2) 
     1412!         ELSE 
     1413!            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * un_adv(i1:i2,j1:j2) 
     1414!         ENDIF 
     1415      ELSE 
     1416         zrhoy = Agrif_Rhoy() 
     1417         ! 
     1418         ubdy(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) &  
     1419                           & / (zrhoy*e2u(i1:i2,j1:j2)) * umask(i1:i2,j1:j2,1) 
     1420         ! 
     1421      ENDIF 
     1422      !  
     1423   END SUBROUTINE interpub2b_const 
     1424 
     1425 
     1426   SUBROUTINE ub2b_cor( ptab, i1, i2, j1, j2, before ) 
     1427      !!---------------------------------------------------------------------- 
     1428      !!                  ***  ROUTINE ub2b_cor  *** 
     1429      !!----------------------------------------------------------------------   
     1430      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     1431      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     1432      LOGICAL                         , INTENT(in   ) ::   before 
     1433      ! 
     1434      INTEGER  :: ji, jj 
     1435      REAL(wp) :: zrhox, zrhoy, zx 
     1436      !!----------------------------------------------------------------------   
     1437      IF( before ) THEN 
     1438         ptab(:,:) = 0._wp 
     1439         DO ji=i1+1,i2-1 
     1440            DO jj=j1+1,j2 
     1441               ptab(ji,jj) = 0.25_wp*( ( vb2_b(ji+1,jj  )*e1v(ji+1,jj  )   &  
     1442                           &            -vb2_b(ji-1,jj  )*e1v(ji-1,jj  ) ) & 
     1443                           &          -( vb2_b(ji+1,jj-1)*e1v(ji+1,jj-1)   & 
     1444                           &            -vb2_b(ji-1,jj-1)*e1v(ji-1,jj-1) ) ) 
     1445            END DO 
     1446         END DO  
     1447      ELSE 
     1448         ! 
     1449         zrhox = Agrif_Rhox()  
     1450         zrhoy = Agrif_Rhoy() 
     1451         DO ji=i1,i2 
     1452            DO jj=j1,j2 
     1453               IF (utint_stage(ji,jj)==0) THEN  
     1454                  zx = 2._wp*MOD(ABS(mig0(ji)-nbghostcells-1), INT(Agrif_Rhox()))/zrhox - 1._wp   
     1455                  ubdy(ji,jj) = ubdy(ji,jj) + 0.25_wp*(1._wp-zx*zx) * ptab(ji,jj) &  
     1456                              &         / zrhoy *r1_e2u(ji,jj) * umask(ji,jj,1)  
     1457                  utint_stage(ji,jj) = 1  
     1458               ENDIF 
     1459            END DO 
     1460         END DO  
     1461         ! 
     1462      ENDIF 
     1463      !  
     1464   END SUBROUTINE ub2b_cor 
     1465 
    12171466 
    12181467   SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before ) 
    … …  
    12521501 
    12531502 
     1503   SUBROUTINE interpvb2b_const( ptab, i1, i2, j1, j2, before ) 
     1504      !!---------------------------------------------------------------------- 
     1505      !!                  ***  ROUTINE interpub2b_const  *** 
     1506      !!----------------------------------------------------------------------   
     1507      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     1508      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     1509      LOGICAL                         , INTENT(in   ) ::   before 
     1510      ! 
     1511      REAL(wp) :: zrhox 
     1512      !!----------------------------------------------------------------------   
     1513      IF( before ) THEN 
     1514!         IF ( ln_bt_fw ) THEN 
     1515            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2) 
     1516!         ELSE 
     1517!            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2) 
     1518!         ENDIF 
     1519      ELSE 
     1520         zrhox = Agrif_Rhox() 
     1521         ! 
     1522         vbdy(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) & 
     1523                           & / (zrhox*e1v(i1:i2,j1:j2)) * vmask(i1:i2,j1:j2,1) 
     1524         ! 
     1525      ENDIF 
     1526      !  
     1527   END SUBROUTINE interpvb2b_const 
     1528 
     1529  
     1530   SUBROUTINE vb2b_cor( ptab, i1, i2, j1, j2, before ) 
     1531      !!---------------------------------------------------------------------- 
     1532      !!                  ***  ROUTINE vb2b_cor  *** 
     1533      !!----------------------------------------------------------------------   
     1534      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     1535      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     1536      LOGICAL                         , INTENT(in   ) ::   before 
     1537      ! 
     1538      INTEGER  :: ji, jj 
     1539      REAL(wp) :: zrhox, zrhoy, zy 
     1540      !!----------------------------------------------------------------------   
     1541      IF( before ) THEN 
     1542         ptab(:,:) = 0._wp 
     1543         DO ji=i1+1,i2-1 
     1544            DO jj=j1+1,j2 
     1545               ptab(ji,jj) = 0.25_wp*( ( ub2_b(ji  ,jj+1)*e2u(ji  ,jj+1)   &  
     1546                           &            -ub2_b(ji  ,jj-1)*e2u(ji  ,jj-1) ) & 
     1547                           &          -( ub2_b(ji-1,jj+1)*e2u(ji-1,jj+1)   & 
     1548                           &            -ub2_b(ji-1,jj-1)*e2u(ji-1,jj-1) ) ) 
     1549            END DO 
     1550         END DO  
     1551      ELSE 
     1552         ! 
     1553         zrhox = Agrif_Rhox()  
     1554         zrhoy = Agrif_Rhoy() 
     1555         DO ji=i1,i2 
     1556            DO jj=j1,j2 
     1557               IF (vtint_stage(ji,jj)==0) THEN  
     1558                  zy = 2._wp*MOD(ABS(mjg0(jj)-nbghostcells-1), INT(Agrif_Rhoy()))/zrhoy - 1._wp   
     1559                  vbdy(ji,jj) = vbdy(ji,jj) + 0.25_wp*(1._wp-zy*zy) * ptab(ji,jj) &  
     1560                              &         / zrhox * r1_e1v(ji,jj) * vmask(ji,jj,1)  
     1561                  vtint_stage(ji,jj) = 1  
     1562               ENDIF 
     1563            END DO 
     1564         END DO  
     1565         ! 
     1566      ENDIF 
     1567      !  
     1568   END SUBROUTINE vb2b_cor 
     1569 
     1570 
    12541571   SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before ) 
    12551572      !!---------------------------------------------------------------------- 
    … …  
    12731590                     WRITE(numout,*) ' Agrif error for e3t_0: parent , child, i, j, k ',  &  
    12741591                     &                 ptab(ji,jj,jk), tmask(ji,jj,jk) * e3t_0(ji,jj,jk), & 
    1275                      &                 mig0(ji), mig0(jj), jk 
     1592                     &                 mig0(ji), mjg0(jj), jk 
    12761593                !     kindic_agr = kindic_agr + 1 
    12771594                  ENDIF 
    … …  
    12831600      !  
    12841601   END SUBROUTINE interpe3t 
     1602 
     1603 
     1604   SUBROUTINE interpe3t0_vremap( ptab, i1, i2, j1, j2, k1, k2, before ) 
     1605      !!---------------------------------------------------------------------- 
     1606      !!                  ***  ROUTINE interpe3t0_vremap  *** 
     1607      !!----------------------------------------------------------------------   
     1608      INTEGER                              , INTENT(in   ) :: i1, i2, j1, j2, k1, k2 
     1609      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab 
     1610      LOGICAL                              , INTENT(in   ) :: before 
     1611      ! 
     1612      INTEGER :: ji, jj, jk 
     1613      REAL(wp) :: zh 
     1614      !!----------------------------------------------------------------------   
     1615      !     
     1616      IF( before ) THEN 
     1617         IF ( ln_zps ) THEN 
     1618            DO jk = k1, k2 
     1619               DO jj = j1, j2 
     1620                  DO ji = i1, i2 
     1621                     ptab(ji, jj, jk) = e3t_1d(jk) 
     1622                  END DO 
     1623               END DO 
     1624            END DO 
     1625         ELSE 
     1626            ptab(i1:i2,j1:j2,k1:k2) = e3t_0(i1:i2,j1:j2,k1:k2) 
     1627         ENDIF 
     1628      ELSE 
     1629         ! 
     1630         DO jk = k1, k2 
     1631            DO jj = j1, j2 
     1632               DO ji = i1, i2 
     1633                  e3t0_parent(ji,jj,jk) = ptab(ji,jj,jk) 
     1634               END DO 
     1635            END DO 
     1636         END DO 
     1637 
     1638         ! Retrieve correct scale factor at the bottom: 
     1639         DO jj = j1, j2 
     1640            DO ji = i1, i2 
     1641               zh = 0._wp 
     1642               DO jk = 1, mbkt_parent(ji, jj)-1 
     1643                  zh = zh + e3t0_parent(ji,jj,jk) 
     1644               END DO 
     1645               e3t0_parent(ji,jj,mbkt_parent(ji,jj)) = ht0_parent(ji, jj) - zh 
     1646            END DO 
     1647         END DO 
     1648          
     1649      ENDIF 
     1650      !  
     1651   END SUBROUTINE interpe3t0_vremap 
     1652 
    12851653 
    12861654   SUBROUTINE interpglamt( ptab, i1, i2, j1, j2, before ) 
    … …  
    14291797   SUBROUTINE interpmbkt( ptab, i1, i2, j1, j2, before ) 
    14301798      !!---------------------------------------------------------------------- 
    1431       !!                  ***  ROUTINE interpsshn  *** 
     1799      !!                  ***  ROUTINE interpmbkt  *** 
    14321800      !!----------------------------------------------------------------------   
    14331801      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
    … …  
    14481816   SUBROUTINE interpht0( ptab, i1, i2, j1, j2, before ) 
    14491817      !!---------------------------------------------------------------------- 
    1450       !!                  ***  ROUTINE interpsshn  *** 
     1818      !!                  ***  ROUTINE interpht0  *** 
    14511819      !!----------------------------------------------------------------------   
    14521820      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
    … …  
    14641832   END SUBROUTINE interpht0 
    14651833 
    1466     
    1467    SUBROUTINE agrif_initts(tabres,i1,i2,j1,j2,k1,k2,m1,m2,before) 
    1468        INTEGER :: i1, i2, j1, j2, k1, k2, m1, m2 
    1469        REAL(wp):: tabres(i1:i2,j1:j2,k1:k2,m1:m2) 
    1470        LOGICAL :: before 
    1471  
    1472        INTEGER :: jm 
    1473  
    1474        IF (before) THEN 
    1475          DO jm=1,jpts 
    1476              tabres(i1:i2,j1:j2,k1:k2,jm) = ts(i1:i2,j1:j2,k1:k2,jm,Kbb_a) 
    1477          END DO 
    1478        ELSE 
    1479          DO jm=1,jpts 
    1480              ts(i1:i2,j1:j2,k1:k2,jm,Kbb_a)=tabres(i1:i2,j1:j2,k1:k2,jm) 
    1481          END DO 
    1482        ENDIF 
    1483    END SUBROUTINE agrif_initts  
    1484  
    1485     
    1486    SUBROUTINE agrif_initssh( ptab, i1, i2, j1, j2, before ) 
    1487       !!---------------------------------------------------------------------- 
    1488       !!                  ***  ROUTINE interpsshn  *** 
    1489       !!----------------------------------------------------------------------   
    1490       INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
    1491       REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
    1492       LOGICAL                         , INTENT(in   ) ::   before 
    1493       ! 
    1494       !!----------------------------------------------------------------------   
    1495       ! 
    1496       IF( before) THEN 
    1497          ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Kbb_a) 
    1498       ELSE 
    1499          ssh(i1:i2,j1:j2,Kbb_a) = ptab(i1:i2,j1:j2)*tmask(i1:i2,j1:j2,1) 
    1500       ENDIF 
    1501       ! 
    1502    END SUBROUTINE agrif_initssh 
     1834   SUBROUTINE Agrif_check_bat( iindic ) 
     1835      !!---------------------------------------------------------------------- 
     1836      !!                  ***  ROUTINE Agrif_check_bat  *** 
     1837      !!----------------------------------------------------------------------   
     1838      INTEGER, INTENT(inout) ::   iindic 
     1839      !! 
     1840      INTEGER :: ji, jj 
     1841      INTEGER  :: istart, iend, jstart, jend, ispon 
     1842      !!----------------------------------------------------------------------   
     1843      ! 
     1844      ! 
     1845      ! --- West --- ! 
     1846      IF(lk_west) THEN 
     1847         ispon  = nn_sponge_len * Agrif_irhox() 
     1848         istart = nn_hls + 2                                  ! halo + land + 1 
     1849         iend   = nn_hls + 1 + nbghostcells + ispon           ! halo + land + nbghostcells + sponge 
     1850         jstart = nn_hls + 2 
     1851         jend   = jpjglo - nn_hls - 1 
     1852         DO ji = mi0(istart), mi1(iend) 
     1853            DO jj = mj0(jstart), mj1(jend) 
     1854               IF ( ABS(ht0_parent(ji,jj)-ht_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1 
     1855            END DO 
     1856            DO jj = mj0(jstart), mj1(jend-1) 
     1857               IF ( ABS(hv0_parent(ji,jj)-hv_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1 
     1858            END DO 
     1859         END DO 
     1860         DO ji = mi0(istart), mi1(iend-1) 
     1861            DO jj = mj0(jstart), mj1(jend) 
     1862               IF ( ABS(hu0_parent(ji,jj)-hu_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1 
     1863            END DO 
     1864         END DO 
     1865      ENDIF 
     1866      ! 
     1867      ! --- East --- ! 
     1868      IF(lk_east) THEN 
     1869         ispon  = nn_sponge_len * Agrif_irhox()  
     1870         istart = jpiglo - ( nn_hls + nbghostcells + ispon )  ! halo + land + nbghostcells + sponge - 1 
     1871         iend   = jpiglo - ( nn_hls + 1 )                     ! halo + land + 1                     - 1 
     1872         jstart = nn_hls + 2 
     1873         jend   = jpjglo - nn_hls - 1  
     1874         DO ji = mi0(istart), mi1(iend) 
     1875            DO jj = mj0(jstart), mj1(jend) 
     1876               IF ( ABS(ht0_parent(ji,jj)-ht_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1 
     1877            END DO 
     1878            DO jj = mj0(jstart), mj1(jend-1) 
     1879               IF ( ABS(hv0_parent(ji,jj)-hv_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1 
     1880            END DO 
     1881         END DO 
     1882         DO ji = mi0(istart+1), mi1(iend-1) 
     1883            DO jj = mj0(jstart), mj1(jend) 
     1884               IF ( ABS(hu0_parent(ji,jj)-hu_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1 
     1885            END DO 
     1886         END DO 
     1887      ENDIF 
     1888      ! 
     1889      ! --- South --- ! 
     1890      IF(lk_south) THEN 
     1891         ispon  = nn_sponge_len * Agrif_irhoy()  
     1892         jstart = nn_hls + 2                                 ! halo + land + 1 
     1893         jend   = nn_hls + 1 + nbghostcells + ispon          ! halo + land + nbghostcells + sponge 
     1894         istart = nn_hls + 2 
     1895         iend   = jpiglo - nn_hls - 1 
     1896         DO jj = mj0(jstart), mj1(jend) 
     1897            DO ji = mi0(istart), mi1(iend) 
     1898               IF ( ABS(ht0_parent(ji,jj)-ht_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1 
     1899            END DO 
     1900            DO ji = mi0(istart), mi1(iend-1) 
     1901               IF ( ABS(hu0_parent(ji,jj)-hu_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1 
     1902            END DO 
     1903         END DO 
     1904         DO jj = mj0(jstart), mj1(jend-1) 
     1905            DO ji = mi0(istart), mi1(iend) 
     1906               IF ( ABS(hv0_parent(ji,jj)-hv_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1 
     1907            END DO 
     1908         END DO 
     1909      ENDIF 
     1910      ! 
     1911      ! --- North --- ! 
     1912      IF(lk_north) THEN 
     1913         ispon  = nn_sponge_len * Agrif_irhoy()  
     1914         jstart = jpjglo - ( nn_hls + nbghostcells + ispon)  ! halo + land + nbghostcells +sponge - 1 
     1915         jend   = jpjglo - ( nn_hls + 1 )                    ! halo + land + 1            - 1 
     1916         istart = nn_hls + 2 
     1917         iend   = jpiglo - nn_hls - 1 
     1918         DO jj = mj0(jstart), mj1(jend) 
     1919            DO ji = mi0(istart), mi1(iend) 
     1920               IF ( ABS(ht0_parent(ji,jj)-ht_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1 
     1921            END DO 
     1922            DO ji = mi0(istart), mi1(iend-1) 
     1923               IF ( ABS(hu0_parent(ji,jj)-hu_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1 
     1924            END DO 
     1925         END DO 
     1926         DO jj = mj0(jstart+1), mj1(jend-1) 
     1927            DO ji = mi0(istart), mi1(iend) 
     1928               IF ( ABS(hv0_parent(ji,jj)-hv_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1 
     1929            END DO 
     1930         END DO 
     1931      ENDIF 
     1932      ! 
     1933   END SUBROUTINE Agrif_check_bat 
    15031934    
    15041935#else 

  • NEMO/trunk/src/NST/agrif_oce_sponge.F90

    r14053 r14086  
    2828   PRIVATE 
    2929 
    30    PUBLIC Agrif_Sponge, Agrif_Sponge_Tra, Agrif_Sponge_Dyn 
     30   PUBLIC Agrif_Sponge, Agrif_Sponge_2d, Agrif_Sponge_Tra, Agrif_Sponge_Dyn 
    3131   PUBLIC interptsn_sponge, interpun_sponge, interpvn_sponge 
     32   PUBLIC interpunb_sponge, interpvnb_sponge 
    3233 
    3334   !! * Substitutions 
    … …  
    4647      !!---------------------------------------------------------------------- 
    4748      REAL(wp) ::   zcoef   ! local scalar 
    48        
     49      INTEGER  :: istart, iend, jstart, jend  
    4950      !!---------------------------------------------------------------------- 
    5051      ! 
    … …  
    5354      zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot()) 
    5455 
    55       CALL Agrif_Sponge 
    5656      Agrif_SpecialValue    = 0._wp 
    5757      Agrif_UseSpecialValue = .TRUE. 
     58      l_vremap              = ln_vert_remap 
    5859      tabspongedone_tsn     = .FALSE. 
    5960      ! 
    60       CALL Agrif_Bc_Variable( tsn_sponge_id, calledweight=zcoef, procname=interptsn_sponge ) 
     61      CALL Agrif_Bc_Variable( ts_sponge_id, calledweight=zcoef, procname=interptsn_sponge ) 
    6162      ! 
    6263      Agrif_UseSpecialValue = .FALSE. 
     64      l_vremap              = .FALSE. 
    6365#endif 
    6466      ! 
    … …  
    8183      Agrif_SpecialValue    = 0._wp 
    8284      Agrif_UseSpecialValue = ln_spc_dyn 
     85      l_vremap              = ln_vert_remap 
    8386      use_sign_north        = .TRUE. 
    8487      sign_north            = -1._wp 
    … …  
    9194      tabspongedone_v = .FALSE. 
    9295      CALL Agrif_Bc_Variable( vn_sponge_id, calledweight=zcoef, procname=interpvn_sponge ) 
     96       
     97      IF ( nn_shift_bar>0 ) THEN ! then split sponge between 2d and 3d 
     98         zcoef = REAL(Agrif_NbStepint(),wp)/REAL(Agrif_rhot()) ! forward tsplit 
     99         tabspongedone_u = .FALSE. 
     100         tabspongedone_v = .FALSE.          
     101         CALL Agrif_Bc_Variable( unb_sponge_id, calledweight=zcoef, procname=interpunb_sponge ) 
     102         ! 
     103         tabspongedone_u = .FALSE. 
     104         tabspongedone_v = .FALSE. 
     105         CALL Agrif_Bc_Variable( vnb_sponge_id, calledweight=zcoef, procname=interpvnb_sponge ) 
     106      ENDIF 
    93107      ! 
    94108      Agrif_UseSpecialValue = .FALSE. 
    95109      use_sign_north        = .FALSE. 
     110      l_vremap              = .FALSE. 
     111      ! 
    96112#endif 
    97113      ! 
    … …  
    110126      REAL(wp) ::   z1_ispongearea, z1_jspongearea 
    111127      REAL(wp), DIMENSION(jpi,jpj) :: ztabramp 
    112 #if defined key_vertical 
    113       REAL(wp), DIMENSION(jpi,jpj) :: ztabrampu 
    114       REAL(wp), DIMENSION(jpi,jpj) :: ztabrampv 
    115 #endif 
    116       REAL(wp), DIMENSION(jpjmax)  :: zmskwest,  zmskeast 
    117       REAL(wp), DIMENSION(jpimax)  :: zmsknorth, zmsksouth 
    118128      !!---------------------------------------------------------------------- 
    119129      ! 
    … …  
    124134      !|            |           |           |           | 
    125135      !fine :     t u t u t u t u t u t u t u t u t u t u t 
    126       !sponge val:0   0   0   1  5/6 4/6 3/6 2/6 1/6  0   0 
     136      !sponge val:0 1   1   1   1  5/6 4/6 3/6 2/6 1/6  0  
    127137      !           |   ghost     | <-- sponge area  -- > | 
    128138      !           |   points    |                       | 
    … …  
    130140 
    131141#if defined SPONGE || defined SPONGE_TOP 
    132       IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN 
    133          ! 
    134          ! Retrieve masks at open boundaries: 
    135  
    136          IF( lk_west ) THEN                             ! --- West --- ! 
    137             ztabramp(:,:) = 0._wp 
    138             ind1 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
    139             DO ji = mi0(ind1), mi1(ind1)                 
    140                ztabramp(ji,:) = ssumask(ji,:) 
    141             END DO 
    142             zmskwest(    1:jpj   ) = MAXVAL(ztabramp(:,:), dim=1) 
    143             zmskwest(jpj+1:jpjmax) = 0._wp 
    144          ENDIF 
    145          IF( lk_east ) THEN                             ! --- East --- ! 
    146             ztabramp(:,:) = 0._wp 
    147             ind1 = jpiglo - ( nn_hls + nbghostcells + 1)   ! halo + land + nbghostcells 
    148             DO ji = mi0(ind1), mi1(ind1)                  
    149                ztabramp(ji,:) = ssumask(ji,:) 
    150             END DO 
    151             zmskeast(    1:jpj   ) = MAXVAL(ztabramp(:,:), dim=1) 
    152             zmskeast(jpj+1:jpjmax) = 0._wp 
    153          ENDIF 
    154          IF( lk_south ) THEN                            ! --- South --- ! 
    155             ztabramp(:,:) = 0._wp 
    156             ind1 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
    157             DO jj = mj0(ind1), mj1(ind1)                  
    158                ztabramp(:,jj) = ssvmask(:,jj) 
    159             END DO 
    160             zmsksouth(    1:jpi   ) = MAXVAL(ztabramp(:,:), dim=2) 
    161             zmsksouth(jpi+1:jpimax) = 0._wp 
    162          ENDIF 
    163          IF( lk_north ) THEN                            ! --- North --- ! 
    164             ztabramp(:,:) = 0._wp 
    165             ind1 = jpjglo - ( nn_hls + nbghostcells + 1)   ! halo + land + nbghostcells 
    166             DO jj = mj0(ind1), mj1(ind1)                  
    167                ztabramp(:,jj) = ssvmask(:,jj) 
    168             END DO 
    169             zmsknorth(    1:jpi   ) = MAXVAL(ztabramp(:,:), dim=2) 
    170             zmsknorth(jpi+1:jpimax) = 0._wp 
    171          ENDIF 
    172  
    173          ! JC: SPONGE MASKING TO BE SORTED OUT: 
    174          zmskwest(:)  = 1._wp 
    175          zmskeast(:)  = 1._wp 
    176          zmsksouth(:) = 1._wp 
    177          zmsknorth(:) = 1._wp 
    178 #if defined key_mpp_mpi 
    179 !         CALL mpp_max( 'AGRIF_sponge', zmskwest(:) , jpjmax ) 
    180 !         CALL mpp_max( 'AGRIF_Sponge', zmskeast(:) , jpjmax ) 
    181 !         CALL mpp_max( 'AGRIF_Sponge', zmsksouth(:), jpimax ) 
    182 !         CALL mpp_max( 'AGRIF_Sponge', zmsknorth(:), jpimax ) 
     142      ! Define ramp from boundaries towards domain interior at F-points 
     143      ! Store it in ztabramp 
     144 
     145      ispongearea    = nn_sponge_len * Agrif_irhox() 
     146      z1_ispongearea = 1._wp / REAL( MAX(ispongearea,1), wp ) 
     147      jspongearea    = nn_sponge_len * Agrif_irhoy() 
     148      z1_jspongearea = 1._wp / REAL( MAX(jspongearea,1), wp ) 
     149          
     150      ztabramp(:,:) = 0._wp 
     151 
     152      IF( lk_west ) THEN                             ! --- West --- ! 
     153         ind1 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
     154         ind2 = nn_hls + 1 + nbghostcells + ispongearea  
     155         DO ji = mi0(ind1), mi1(ind2)    
     156            DO jj = 1, jpj                
     157               ztabramp(ji,jj) =                       REAL(ind2 - mig(ji), wp) * z1_ispongearea 
     158            END DO 
     159         END DO 
     160         ! ghost cells: 
     161         ind1 = 1 
     162         ind2 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
     163         DO ji = mi0(ind1), mi1(ind2)    
     164            DO jj = 1, jpj                
     165               ztabramp(ji,jj) = 1._wp 
     166            END DO 
     167         END DO 
     168      ENDIF 
     169      IF( lk_east ) THEN                             ! --- East --- ! 
     170         ind1 = jpiglo - ( nn_hls + nbghostcells ) - ispongearea - 1 
     171         ind2 = jpiglo - ( nn_hls + nbghostcells ) - 1    ! halo + land + nbghostcells - 1 
     172         DO ji = mi0(ind1), mi1(ind2) 
     173            DO jj = 1, jpj 
     174               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mig(ji) - ind1, wp) * z1_ispongearea )  
     175            END DO 
     176         END DO 
     177         ! ghost cells: 
     178         ind1 = jpiglo - ( nn_hls + nbghostcells ) - 1    ! halo + land + nbghostcells - 1 
     179         ind2 = jpiglo - 1 
     180         DO ji = mi0(ind1), mi1(ind2) 
     181            DO jj = 1, jpj 
     182               ztabramp(ji,jj) = 1._wp 
     183            END DO 
     184         END DO 
     185      ENDIF       
     186      IF( lk_south ) THEN                            ! --- South --- ! 
     187         ind1 = nn_hls + 1 + nbghostcells                 ! halo + land + nbghostcells 
     188         ind2 = nn_hls + 1 + nbghostcells + jspongearea  
     189         DO jj = mj0(ind1), mj1(ind2)  
     190            DO ji = 1, jpi 
     191               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(ind2 - mjg(jj), wp) * z1_jspongearea ) 
     192            END DO 
     193         END DO 
     194         ! ghost cells: 
     195         ind1 = 1 
     196         ind2 = nn_hls + 1 + nbghostcells                 ! halo + land + nbghostcells 
     197         DO jj = mj0(ind1), mj1(ind2)  
     198            DO ji = 1, jpi 
     199               ztabramp(ji,jj) = 1._wp 
     200            END DO 
     201         END DO 
     202      ENDIF 
     203      IF( lk_north ) THEN                            ! --- North --- ! 
     204         ind1 = jpjglo - ( nn_hls + nbghostcells ) - jspongearea - 1 
     205         ind2 = jpjglo - ( nn_hls + nbghostcells ) - 1    ! halo + land + nbghostcells - 1 
     206         DO jj = mj0(ind1), mj1(ind2) 
     207            DO ji = 1, jpi 
     208               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mjg(jj) - ind1, wp) * z1_jspongearea )  
     209            END DO 
     210         END DO 
     211         ! ghost cells: 
     212         ind1 = jpjglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
     213         ind2 = jpjglo 
     214         DO jj = mj0(ind1), mj1(ind2) 
     215            DO ji = 1, jpi 
     216               ztabramp(ji,jj) = 1._wp 
     217            END DO 
     218         END DO 
     219      ENDIF       
     220      ! 
     221      ! Tracers 
     222      fspu(:,:) = 0._wp 
     223      fspv(:,:) = 0._wp 
     224      DO_2D( 0, 0, 0, 0 ) 
     225         fspu(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji,jj-1) ) * ssumask(ji,jj) 
     226         fspv(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji-1,jj) ) * ssvmask(ji,jj) 
     227      END_2D 
     228 
     229      ! Dynamics 
     230      fspt(:,:) = 0._wp 
     231      fspf(:,:) = 0._wp 
     232      DO_2D( 0, 0, 0, 0 ) 
     233         fspt(ji,jj) = 0.25_wp * ( ztabramp(ji  ,jj  ) + ztabramp(ji-1,jj  ) & 
     234                               &  +ztabramp(ji  ,jj-1) + ztabramp(ji-1,jj-1) ) * ssmask(ji,jj) 
     235         fspf(ji,jj) = ztabramp(ji,jj) * ssvmask(ji,jj) * ssvmask(ji,jj+1) 
     236      END_2D 
     237       
     238      CALL lbc_lnk_multi( 'agrif_Sponge', fspu, 'U', 1._wp, fspv, 'V', 1._wp, fspt, 'T', 1._wp, fspf, 'F', 1._wp ) 
     239      ! 
     240      ! Remove vertical interpolation where not needed: 
     241      ! (A null value in mbkx arrays does the job) 
     242      WHERE (fspu(:,:) == 0._wp) mbku_parent(:,:) = 0 
     243      WHERE (fspv(:,:) == 0._wp) mbkv_parent(:,:) = 0 
     244      WHERE (fspt(:,:) == 0._wp) mbkt_parent(:,:) = 0 
     245      ! 
    183246#endif 
    184  
    185          ! Define ramp from boundaries towards domain interior at T-points 
    186          ! Store it in ztabramp 
    187  
    188          ispongearea    = nn_sponge_len * Agrif_irhox() 
    189          z1_ispongearea = 1._wp / REAL( ispongearea, wp ) 
    190          jspongearea    = nn_sponge_len * Agrif_irhoy() 
    191          z1_jspongearea = 1._wp / REAL( jspongearea, wp ) 
    192           
    193          ztabramp(:,:) = 0._wp 
    194  
    195          ! Trick to remove sponge in 2DV domains: 
    196          IF ( nbcellsx <= 3 ) ispongearea = -1 
    197          IF ( nbcellsy <= 3 ) jspongearea = -1 
    198  
    199          IF( lk_west ) THEN                             ! --- West --- ! 
    200             ind1 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
    201             ind2 = nn_hls + 1 + nbghostcells + ispongearea  
    202             DO ji = mi0(ind1), mi1(ind2)    
    203                DO jj = 1, jpj                
    204                   ztabramp(ji,jj) =                       REAL(ind2 - mig(ji), wp) * z1_ispongearea   * zmskwest(jj) 
    205                END DO 
    206             END DO 
    207             ! ghost cells: 
    208             ind1 = 1 
    209             ind2 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
    210             DO ji = mi0(ind1), mi1(ind2)    
    211                DO jj = 1, jpj                
    212                   ztabramp(ji,jj) = zmskwest(jj) 
    213                END DO 
    214             END DO 
    215          ENDIF 
    216          IF( lk_east ) THEN                             ! --- East --- ! 
    217             ind1 = jpiglo - ( nn_hls + nbghostcells ) - ispongearea 
    218             ind2 = jpiglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    219             DO ji = mi0(ind1), mi1(ind2) 
    220                DO jj = 1, jpj 
    221                   ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mig(ji) - ind1, wp) * z1_ispongearea ) * zmskeast(jj) 
    222                END DO 
    223             END DO 
    224             ! ghost cells: 
    225             ind1 = jpiglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    226             ind2 = jpiglo 
    227             DO ji = mi0(ind1), mi1(ind2) 
    228                DO jj = 1, jpj 
    229                   ztabramp(ji,jj) = zmskeast(jj) 
    230                END DO 
    231             END DO 
    232          ENDIF       
    233          IF( lk_south ) THEN                            ! --- South --- ! 
    234             ind1 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
    235             ind2 = nn_hls + 1 + nbghostcells + jspongearea  
    236             DO jj = mj0(ind1), mj1(ind2)  
    237                DO ji = 1, jpi 
    238                   ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(ind2 - mjg(jj), wp) * z1_jspongearea ) * zmsksouth(ji) 
    239                END DO 
    240             END DO 
    241             ! ghost cells: 
    242             ind1 = 1 
    243             ind2 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
    244             DO jj = mj0(ind1), mj1(ind2)  
    245                DO ji = 1, jpi 
    246                   ztabramp(ji,jj) = zmsksouth(ji) 
    247                END DO 
    248             END DO 
    249          ENDIF 
    250          IF( lk_north ) THEN                            ! --- North --- ! 
    251             ind1 = jpjglo - ( nn_hls + nbghostcells ) - jspongearea 
    252             ind2 = jpjglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    253             DO jj = mj0(ind1), mj1(ind2) 
    254                DO ji = 1, jpi 
    255                   ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mjg(jj) - ind1, wp) * z1_jspongearea ) * zmsknorth(ji) 
    256                END DO 
    257             END DO 
    258             ! ghost cells: 
    259             ind1 = jpjglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    260             ind2 = jpjglo 
    261             DO jj = mj0(ind1), mj1(ind2) 
    262                DO ji = 1, jpi 
    263                   ztabramp(ji,jj) = zmsknorth(ji) 
    264                END DO 
    265             END DO 
    266          ENDIF 
    267          ! 
     247      ! 
     248   END SUBROUTINE Agrif_Sponge 
     249 
     250 
     251   SUBROUTINE Agrif_Sponge_2d 
     252      !!---------------------------------------------------------------------- 
     253      !!                 *** ROUTINE  Agrif_Sponge_2d *** 
     254      !!---------------------------------------------------------------------- 
     255      INTEGER  ::   ji, jj, ind1, ind2, ishift, jshift 
     256      INTEGER  ::   ispongearea, jspongearea 
     257      REAL(wp) ::   z1_ispongearea, z1_jspongearea 
     258      REAL(wp), DIMENSION(jpi,jpj) :: ztabramp 
     259      !!---------------------------------------------------------------------- 
     260      ! 
     261      ! Sponge 1d example with: 
     262      !      iraf = 3 ; nbghost = 3 ; nn_sponge_len = 2 
     263      !                         
     264      !coarse :     U     T     U     T     U     T     U 
     265      !|            |           |           |           | 
     266      !fine :     t u t u t u t u t u t u t u t u t u t u t 
     267      !sponge val:0 1   1   1   1  5/6 4/6 3/6 2/6 1/6  0  
     268      !           |   ghost     | <-- sponge area  -- > | 
     269      !           |   points    |                       | 
     270      !                         |--> dynamical interface 
     271 
     272#if defined SPONGE || defined SPONGE_TOP 
     273      ! Define ramp from boundaries towards domain interior at F-points 
     274      ! Store it in ztabramp 
     275 
     276      ispongearea    = nn_sponge_len * Agrif_irhox() 
     277      z1_ispongearea = 1._wp / REAL( MAX(ispongearea,1), wp ) 
     278      jspongearea    = nn_sponge_len * Agrif_irhoy() 
     279      z1_jspongearea = 1._wp / REAL( MAX(jspongearea,1), wp ) 
     280      ishift = nn_shift_bar * Agrif_irhox() 
     281      jshift = nn_shift_bar * Agrif_irhoy() 
     282         
     283      ztabramp(:,:) = 0._wp 
     284 
     285      IF( lk_west ) THEN                             ! --- West --- ! 
     286         ind1 = nn_hls + 1 + nbghostcells + ishift 
     287         ind2 = nn_hls + 1 + nbghostcells + ishift + ispongearea  
     288         DO ji = mi0(ind1), mi1(ind2)    
     289            DO jj = 1, jpj                
     290               ztabramp(ji,jj) =                       REAL(ind2 - mig(ji), wp) * z1_ispongearea 
     291            END DO 
     292         END DO 
     293         ! ghost cells: 
     294         ind1 = 1 
     295         ind2 = nn_hls + 1 + nbghostcells + ishift               ! halo + land + nbghostcells 
     296         DO ji = mi0(ind1), mi1(ind2)    
     297            DO jj = 1, jpj                
     298               ztabramp(ji,jj) = 1._wp 
     299            END DO 
     300         END DO 
    268301      ENDIF 
    269  
     302      IF( lk_east ) THEN                             ! --- East --- ! 
     303         ind1 = jpiglo - ( nn_hls + nbghostcells + ishift) - ispongearea - 1 
     304         ind2 = jpiglo - ( nn_hls + nbghostcells + ishift) - 1    ! halo + land + nbghostcells - 1 
     305         DO ji = mi0(ind1), mi1(ind2) 
     306            DO jj = 1, jpj 
     307               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mig(ji) - ind1, wp) * z1_ispongearea )  
     308            END DO 
     309         END DO 
     310         ! ghost cells: 
     311         ind1 = jpiglo - ( nn_hls + nbghostcells + ishift) - 1    ! halo + land + nbghostcells - 1 
     312         ind2 = jpiglo - 1 
     313         DO ji = mi0(ind1), mi1(ind2) 
     314            DO jj = 1, jpj 
     315               ztabramp(ji,jj) = 1._wp 
     316            END DO 
     317         END DO 
     318      ENDIF       
     319      IF( lk_south ) THEN                            ! --- South --- ! 
     320         ind1 = nn_hls + 1 + nbghostcells + jshift                ! halo + land + nbghostcells 
     321         ind2 = nn_hls + 1 + nbghostcells + jshift + jspongearea  
     322         DO jj = mj0(ind1), mj1(ind2)  
     323            DO ji = 1, jpi 
     324               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(ind2 - mjg(jj), wp) * z1_jspongearea ) 
     325            END DO 
     326         END DO 
     327         ! ghost cells: 
     328         ind1 = 1 
     329         ind2 = nn_hls + 1 + nbghostcells + jshift                ! halo + land + nbghostcells 
     330         DO jj = mj0(ind1), mj1(ind2)  
     331            DO ji = 1, jpi 
     332               ztabramp(ji,jj) = 1._wp 
     333            END DO 
     334         END DO 
     335      ENDIF 
     336      IF( lk_north ) THEN                            ! --- North --- ! 
     337         ind1 = jpjglo - ( nn_hls + nbghostcells + jshift) - jspongearea - 1 
     338         ind2 = jpjglo - ( nn_hls + nbghostcells + jshift) - 1    ! halo + land + nbghostcells - 1 
     339         DO jj = mj0(ind1), mj1(ind2) 
     340            DO ji = 1, jpi 
     341               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mjg(jj) - ind1, wp) * z1_jspongearea )  
     342            END DO 
     343         END DO 
     344         ! ghost cells: 
     345         ind1 = jpjglo - ( nn_hls + nbghostcells + jshift)      ! halo + land + nbghostcells - 1 
     346         ind2 = jpjglo 
     347         DO jj = mj0(ind1), mj1(ind2) 
     348            DO ji = 1, jpi 
     349               ztabramp(ji,jj) = 1._wp 
     350            END DO 
     351         END DO 
     352      ENDIF 
     353      ! 
    270354      ! Tracers 
    271       IF( .NOT. spongedoneT ) THEN 
    272          fspu(:,:) = 0._wp 
    273          fspv(:,:) = 0._wp 
    274          DO_2D( 0, 0, 0, 0 ) 
    275             fspu(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji+1,jj  ) ) * ssumask(ji,jj) 
    276             fspv(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji  ,jj+1) ) * ssvmask(ji,jj) 
     355      fspu_2d(:,:) = 0._wp 
     356      fspv_2d(:,:) = 0._wp 
     357      DO_2D( 0, 0, 0, 0 ) 
     358         fspu_2d(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji,jj-1) ) * ssumask(ji,jj) 
     359         fspv_2d(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji-1,jj) ) * ssvmask(ji,jj) 
     360      END_2D 
     361 
     362      ! Dynamics 
     363      fspt_2d(:,:) = 0._wp 
     364      fspf_2d(:,:) = 0._wp 
     365      DO_2D( 0, 0, 0, 0 ) 
     366         fspt_2d(ji,jj) = 0.25_wp * ( ztabramp(ji  ,jj  ) + ztabramp(ji-1,jj  ) & 
     367                               &  +ztabramp(ji  ,jj-1) + ztabramp(ji-1,jj-1) ) * ssmask(ji,jj) 
     368         fspf_2d(ji,jj) = ztabramp(ji,jj) * ssvmask(ji,jj) * ssvmask(ji,jj+1) 
    277369         END_2D 
    278       ENDIF 
    279  
    280       ! Dynamics 
    281       IF( .NOT. spongedoneU ) THEN 
    282          fspt(:,:) = 0._wp 
    283          fspf(:,:) = 0._wp 
    284          DO_2D( 0, 0, 0, 0 ) 
    285             fspt(ji,jj) = ztabramp(ji,jj) * ssmask(ji,jj) 
    286             fspf(ji,jj) = 0.25_wp * ( ztabramp(ji  ,jj  ) + ztabramp(ji  ,jj+1)   & 
    287                                   &  +ztabramp(ji+1,jj+1) + ztabramp(ji+1,jj  ) ) & 
    288                                   &  * ssvmask(ji,jj) * ssvmask(ji,jj+1) 
    289          END_2D 
    290       ENDIF 
    291        
    292       IF( .NOT. spongedoneT .AND. .NOT. spongedoneU ) THEN 
    293          CALL lbc_lnk_multi( 'agrif_Sponge', fspu, 'U', 1._wp, fspv, 'V', 1._wp, fspt, 'T', 1._wp, fspf, 'F', 1._wp ) 
    294          spongedoneT = .TRUE. 
    295          spongedoneU = .TRUE. 
    296       ENDIF 
    297       IF( .NOT. spongedoneT ) THEN 
    298          CALL lbc_lnk_multi( 'agrif_Sponge', fspu, 'U', 1._wp, fspv, 'V', 1._wp ) 
    299          spongedoneT = .TRUE. 
    300       ENDIF 
    301       IF( .NOT. spongedoneT .AND. .NOT. spongedoneU ) THEN 
    302          CALL lbc_lnk_multi( 'agrif_Sponge', fspt, 'T', 1._wp, fspf, 'F', 1._wp ) 
    303          spongedoneU = .TRUE. 
    304       ENDIF 
    305  
    306 #if defined key_vertical 
    307       ! Remove vertical interpolation where not needed: 
    308       DO_2D( 0, 0, 0, 0 ) 
    309          IF ((fspu(ji-1,jj)==0._wp).AND.(fspu(ji,jj)==0._wp).AND. & 
    310          &   (fspv(ji,jj-1)==0._wp).AND.(fspv(ji,jj)==0._wp)) mbkt_parent(ji,jj) = 0 
    311 ! 
    312          IF ((fspt(ji+1,jj)==0._wp).AND.(fspt(ji,jj)==0._wp).AND. & 
    313          &   (fspf(ji,jj-1)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbku_parent(ji,jj) = 0 
    314 ! 
    315          IF ((fspt(ji,jj+1)==0._wp).AND.(fspt(ji,jj)==0._wp).AND. & 
    316          &   (fspf(ji-1,jj)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbkv_parent(ji,jj) = 0 
    317 ! 
    318          IF ( ssmask(ji,jj) == 0._wp) mbkt_parent(ji,jj) = 0 
    319          IF (ssumask(ji,jj) == 0._wp) mbku_parent(ji,jj) = 0 
    320          IF (ssvmask(ji,jj) == 0._wp) mbkv_parent(ji,jj) = 0 
    321       END_2D 
    322       ! 
    323       ztabramp (:,:) = REAL( mbkt_parent(:,:), wp ) 
    324       ztabrampu(:,:) = REAL( mbku_parent(:,:), wp ) 
    325       ztabrampv(:,:) = REAL( mbkv_parent(:,:), wp ) 
    326       CALL lbc_lnk_multi( 'Agrif_Sponge', ztabramp, 'T', 1._wp, ztabrampu, 'U', 1._wp, ztabrampv, 'V', 1._wp ) 
    327       mbkt_parent(:,:) = NINT( ztabramp (:,:) ) 
    328       mbku_parent(:,:) = NINT( ztabrampu(:,:) ) 
    329       mbkv_parent(:,:) = NINT( ztabrampv(:,:) ) 
     370      CALL lbc_lnk_multi( 'agrif_Sponge_2d', fspu_2d, 'U', 1._wp, fspv_2d, 'V', 1._wp, fspt_2d, 'T', 1._wp, fspf_2d, 'F', 1._wp ) 
     371      ! 
    330372#endif 
    331373      ! 
    332 #endif 
    333       ! 
    334    END SUBROUTINE Agrif_Sponge 
    335  
    336     
    337    SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
     374   END SUBROUTINE Agrif_Sponge_2d 
     375 
     376 
     377   SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before)  
    338378      !!---------------------------------------------------------------------- 
    339379      !!                 *** ROUTINE interptsn_sponge *** 
    … …  
    346386      INTEGER  ::   iku, ikv 
    347387      REAL(wp) :: ztsa, zabe1, zabe2, zbtr, zhtot 
    348       REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv 
     388      REAL(wp), DIMENSION(i1-1:i2,j1-1:j2,jpk) :: ztu, ztv 
    349389      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff 
    350390      ! vertical interpolation: 
    351391      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child 
    352       REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin 
    353       REAL(wp), DIMENSION(k1:k2) :: h_in 
    354       REAL(wp), DIMENSION(1:jpk) :: h_out 
     392      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin, tabin_i 
     393      REAL(wp), DIMENSION(k1:k2) :: z_in, z_in_i, h_in_i 
     394      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out 
    355395      INTEGER :: N_in, N_out 
    356396      !!---------------------------------------------------------------------- 
    … …  
    367407         END DO 
    368408 
    369 # if defined key_vertical 
    370         ! Interpolate thicknesses 
    371         ! Warning: these are masked, hence extrapolated prior interpolation. 
    372         DO jk=k1,k2 
    373            DO jj=j1,j2 
    374               DO ji=i1,i2 
    375                   tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kbb_a) 
    376               END DO 
    377            END DO 
    378         END DO 
    379  
    380         ! Extrapolate thicknesses in partial bottom cells: 
    381         ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
    382         IF (ln_zps) THEN 
    383            DO jj=j1,j2 
    384               DO ji=i1,i2 
    385                   jk = mbkt(ji,jj) 
    386                   tabres(ji,jj,jk,jpts+1) = 0._wp 
    387               END DO 
    388            END DO            
    389         END IF 
    390       
    391         ! Save ssh at last level: 
    392         IF (.NOT.ln_linssh) THEN 
    393            tabres(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kbb_a)*tmask(i1:i2,j1:j2,1)  
    394         ELSE 
    395            tabres(i1:i2,j1:j2,k2,jpts+1) = 0._wp 
    396         END IF       
    397 # endif 
    398  
    399       ELSE    
    400          ! 
    401 # if defined key_vertical 
    402  
    403          IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp 
    404  
    405          DO jj=j1,j2 
    406             DO ji=i1,i2 
    407                tabres_child(ji,jj,:,:) = 0._wp  
    408                N_in = mbkt_parent(ji,jj) 
    409                zhtot = 0._wp 
    410                DO jk=1,N_in !k2 = jpk of parent grid 
    411                   IF (jk==N_in) THEN 
    412                      h_in(jk) = ht0_parent(ji,jj) + tabres(ji,jj,k2,n2) - zhtot 
     409         IF ( l_vremap.OR.ln_zps ) THEN 
     410 
     411            ! Fill cell depths (i.e. gdept) to be interpolated 
     412            ! Warning: these are masked, hence extrapolated prior interpolation. 
     413            DO jj=j1,j2 
     414               DO ji=i1,i2 
     415                  tabres(ji,jj,k1,jpts+1) = 0.5_wp * tmask(ji,jj,k1) * e3t(ji,jj,k1,Kbb_a) 
     416                  DO jk=k1+1,k2 
     417                     tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * & 
     418                        & ( tabres(ji,jj,jk-1,jpts+1) + 0.5_wp * (e3t(ji,jj,jk-1,Kbb_a)+e3t(ji,jj,jk,Kbb_a)) ) 
     419                  END DO 
     420               END DO 
     421            END DO 
     422 
     423            ! Save ssh at last level: 
     424            IF ( .NOT.ln_linssh ) THEN 
     425               tabres(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kbb_a)*tmask(i1:i2,j1:j2,1)  
     426            END IF   
     427     
     428         END IF 
     429 
     430      ELSE  
     431         ! 
     432         IF ( l_vremap ) THEN 
     433 
     434            IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp 
     435 
     436            DO jj=j1,j2 
     437               DO ji=i1,i2 
     438 
     439                  tabres_child(ji,jj,:,:) = 0._wp  
     440                  ! Build vertical grids: 
     441                  N_in = mbkt_parent(ji,jj) 
     442                  ! Input grid (account for partial cells if any): 
     443                  DO jk=1,N_in 
     444                     z_in(jk) = tabres(ji,jj,jk,n2) - tabres(ji,jj,k2,n2) 
     445                     tabin(jk,1:jpts) = tabres(ji,jj,jk,1:jpts) 
     446                  END DO 
     447                   
     448                  ! Intermediate grid: 
     449                  DO jk = 1, N_in 
     450                     h_in_i(jk) = e3t0_parent(ji,jj,jk) * &  
     451                       &       (1._wp + tabres(ji,jj,k2,n2)/(ht0_parent(ji,jj)*ssmask(ji,jj) + 1._wp - ssmask(ji,jj))) 
     452                  END DO 
     453                  z_in_i(1) = 0.5_wp * h_in_i(1) 
     454                  DO jk=2,N_in 
     455                     z_in_i(jk) = z_in_i(jk-1) + 0.5_wp * ( h_in_i(jk) + h_in_i(jk-1) ) 
     456                  END DO 
     457                  z_in_i(1:N_in) = z_in_i(1:N_in)  - tabres(ji,jj,k2,n2) 
     458 
     459                  ! Output (Child) grid: 
     460                  N_out = mbkt(ji,jj) 
     461                  DO jk=1,N_out 
     462                     h_out(jk) = e3t(ji,jj,jk,Kbb_a) 
     463                  END DO 
     464                  z_out(1) = 0.5_wp * h_out(1) 
     465                  DO jk=2,N_out 
     466                     z_out(jk) = z_out(jk-1) + 0.5_wp * ( h_out(jk)+h_out(jk-1) ) 
     467                  END DO 
     468                  IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out)  - ssh(ji,jj,Kbb_a) 
     469 
     470                  ! Account for small differences in the free-surface 
     471                  IF ( sum(h_out(1:N_out)) > sum(h_in_i(1:N_in) )) THEN 
     472                     h_out(1) = h_out(1)  - ( sum(h_out(1:N_out))-sum(h_in_i(1:N_in)) ) 
    413473                  ELSE 
    414                      h_in(jk) = tabres(ji,jj,jk,n2) 
     474                     h_in_i(1)= h_in_i(1) - ( sum(h_in_i(1:N_in))-sum(h_out(1:N_out)) ) 
     475                  END IF 
     476                  IF (N_in*N_out > 0) THEN 
     477                     CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),tabin_i(1:N_in,1:jpts),z_in_i(1:N_in),N_in,N_in,jpts) 
     478                     CALL reconstructandremap(tabin_i(1:N_in,1:jpts),h_in_i(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts) 
     479!                     CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),z_out(1:N_in),N_in,N_out,jpts)   
    415480                  ENDIF 
    416                   zhtot = zhtot + h_in(jk) 
    417                   tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1) 
    418                END DO 
    419                N_out = 0 
    420                DO jk=1,jpk ! jpk of child grid 
    421                   IF (tmask(ji,jj,jk) == 0) EXIT  
    422                   N_out = N_out + 1 
    423                   h_out(jk) = e3t(ji,jj,jk,Kbb_a) !Child grid scale factors. Could multiply by e1e2t here instead of division above 
    424                END DO 
    425  
    426                ! Account for small differences in free-surface 
    427                IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 
    428                   h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 
    429                ELSE 
    430                   h_in(1)   = h_in(1)   - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 
    431                ENDIF 
    432                IF (N_in*N_out > 0) THEN 
    433                   CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts) 
    434                ENDIF 
    435             END DO 
    436          END DO 
    437 # endif 
    438  
    439          DO jj=j1,j2 
    440             DO ji=i1,i2 
    441                DO jk=1,jpkm1 
    442 # if defined key_vertical 
    443                   tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres_child(ji,jj,jk,1:jpts)) * tmask(ji,jj,jk) 
    444 # else 
    445                   tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres(ji,jj,jk,1:jpts))*tmask(ji,jj,jk) 
    446 # endif 
    447                END DO 
    448             END DO 
    449          END DO 
     481               END DO 
     482            END DO 
     483 
     484            DO jj=j1,j2 
     485               DO ji=i1,i2 
     486                  DO jk=1,jpkm1 
     487                     tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres_child(ji,jj,jk,1:jpts)) * tmask(ji,jj,jk) 
     488                  END DO 
     489               END DO 
     490            END DO 
     491 
     492         ELSE 
     493 
     494            IF ( Agrif_Parent(ln_zps) ) THEN ! Account for partial cells 
     495 
     496               DO jj=j1,j2 
     497                  DO ji=i1,i2 
     498                     ! 
     499                     N_in  = mbkt(ji,jj)  
     500                     N_out = mbkt(ji,jj)  
     501                     z_in(1) = tabres(ji,jj,1,n2) 
     502                     tabin(1,1:jpts) = tabres(ji,jj,1,1:jpts) 
     503                     DO jk=2, N_in 
     504                        z_in(jk) = tabres(ji,jj,jk,n2) 
     505                        tabin(jk,1:jpts) = tabres(ji,jj,jk,1:jpts) 
     506                     END DO  
     507                     IF (.NOT.ln_linssh) z_in(1:N_in) = z_in(1:N_in) - tabres(ji,jj,k2,n2) 
     508 
     509                     z_out(1) = 0.5_wp * e3t(ji,jj,1,Kbb_a) 
     510                     DO jk=2, N_out 
     511                        z_out(jk) = z_out(jk-1) + 0.5_wp * (e3t(ji,jj,jk-1,Kbb_a) + e3t(ji,jj,jk,Kbb_a))  
     512                     END DO  
     513                     IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out) - ssh(ji,jj,Kbb_a) 
     514 
     515                     CALL remap_linear(tabin(1:N_in,1:jpts), z_in(1:N_in), tabres(ji,jj,1:N_out,1:jpts), & 
     516                                         &   z_out(1:N_out), N_in, N_out, jpts) 
     517                  END DO 
     518               END DO 
     519            ENDIF 
     520 
     521            DO jj=j1,j2 
     522               DO ji=i1,i2 
     523                  DO jk=1,jpkm1 
     524                     tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres(ji,jj,jk,1:jpts))*tmask(ji,jj,jk) 
     525                  END DO 
     526               END DO 
     527            END DO 
     528 
     529         END IF 
    450530 
    451531         DO jn = 1, jpts             
    452532            DO jk = 1, jpkm1 
    453                ztu(i1:i2,j1:j2,jk) = 0._wp 
     533               ztu(i1-1:i2,j1-1:j2,jk) = 0._wp 
    454534               DO jj = j1,j2 
    455535                  DO ji = i1,i2-1 
    456                      zabe1 = rn_sponge_tra * r1_Dt * fspu(ji,jj) * umask(ji,jj,jk) * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) 
    457                      ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) )  
    458                   END DO 
    459                END DO 
    460                ztv(i1:i2,j1:j2,jk) = 0._wp 
     536                     zabe1 = rn_sponge_tra * r1_Dt * umask(ji,jj,jk) * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) 
     537                     ztu(ji,jj,jk) = zabe1 * fspu(ji,jj) * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) )  
     538                  END DO 
     539               END DO 
     540               ztv(i1-1:i2,j1-1:j2,jk) = 0._wp 
    461541               DO ji = i1,i2 
    462542                  DO jj = j1,j2-1 
    463                      zabe2 = rn_sponge_tra * r1_Dt * fspv(ji,jj) * vmask(ji,jj,jk) * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a) 
    464                      ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
     543                     zabe2 = rn_sponge_tra * r1_Dt * vmask(ji,jj,jk) * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a) 
     544                     ztv(ji,jj,jk) = zabe2 * fspv(ji,jj) * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
    465545                  END DO 
    466546               END DO 
    … …  
    479559            END DO 
    480560            ! 
     561! JC: there is something wrong with the Laplacian in corners 
    481562            DO jk = 1, jpkm1 
    482                DO jj = j1+1,j2-1 
    483                   DO ji = i1+1,i2-1 
     563               DO jj = j1,j2 
     564                  DO ji = i1,i2 
    484565                     IF (.NOT. tabspongedone_tsn(ji,jj)) THEN  
    485566                        zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a) 
    486567                        ! horizontal diffusive trends 
    487                         ztsa = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) & 
    488                              &  - rn_trelax_tra * r1_Dt * fspt(ji,jj) * tsbdiff(ji,jj,jk,jn)  
     568                        ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   &  
     569                          &           + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) & 
     570                          &   - rn_trelax_tra * r1_Dt * fspt(ji,jj) * tsbdiff(ji,jj,jk,jn) 
    489571                        ! add it to the general tracer trends 
    490572                        ts(ji,jj,jk,jn,Krhs_a) = ts(ji,jj,jk,jn,Krhs_a) + ztsa 
    … …  
    492574                  END DO 
    493575               END DO 
     576 
    494577            END DO 
    495578            ! 
    496579         END DO 
    497580         ! 
    498          tabspongedone_tsn(i1+1:i2-1,j1+1:j2-1) = .TRUE. 
     581         tabspongedone_tsn(i1:i2,j1:j2) = .TRUE. 
    499582         ! 
    500583      ENDIF 
    … …  
    529612               DO ji=i1,i2 
    530613                  tabres(ji,jj,jk,m1) = uu(ji,jj,jk,Kbb_a) 
    531 # if defined key_vertical 
    532                   tabres(ji,jj,jk,m2) = e3u(ji,jj,jk,Kbb_a)*umask(ji,jj,jk) 
    533 # endif 
    534                END DO 
    535             END DO 
    536          END DO 
    537  
    538 # if defined key_vertical 
    539          ! Extrapolate thicknesses in partial bottom cells: 
    540          ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
    541          IF (ln_zps) THEN 
     614               END DO 
     615            END DO 
     616         END DO 
     617 
     618         IF ( l_vremap ) THEN 
     619 
     620            DO jk=k1,k2 
     621               DO jj=j1,j2 
     622                  DO ji=i1,i2 
     623                     tabres(ji,jj,jk,m2) = e3u(ji,jj,jk,Kbb_a)*umask(ji,jj,jk) 
     624                  END DO 
     625               END DO 
     626            END DO 
     627 
     628            ! Extrapolate thicknesses in partial bottom cells: 
     629            ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
     630            IF (ln_zps) THEN 
     631               DO jj=j1,j2 
     632                  DO ji=i1,i2 
     633                     jk = mbku(ji,jj) 
     634                     tabres(ji,jj,jk,m2) = 0._wp 
     635                  END DO 
     636               END DO            
     637            END IF 
     638            ! Save ssh at last level: 
     639            tabres(i1:i2,j1:j2,k2,m2) = 0._wp 
     640            IF (.NOT.ln_linssh) THEN 
     641               ! This vertical sum below should be replaced by the sea-level at U-points (optimization): 
     642               DO jk=1,jpk 
     643                  tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3u(i1:i2,j1:j2,jk,Kbb_a) * umask(i1:i2,j1:j2,jk) 
     644               END DO 
     645               tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hu_0(i1:i2,j1:j2) 
     646            END IF  
     647         END IF 
     648 
     649      ELSE 
     650 
     651         IF ( l_vremap ) THEN 
     652 
     653            IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp 
     654 
    542655            DO jj=j1,j2 
    543656               DO ji=i1,i2 
    544                   jk = mbku(ji,jj) 
    545                   tabres(ji,jj,jk,m2) = 0._wp 
    546                END DO 
    547             END DO            
    548          END IF 
    549         ! Save ssh at last level: 
    550         tabres(i1:i2,j1:j2,k2,m2) = 0._wp 
    551         IF (.NOT.ln_linssh) THEN 
    552            ! This vertical sum below should be replaced by the sea-level at U-points (optimization): 
    553            DO jk=1,jpk 
    554               tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3u(i1:i2,j1:j2,jk,Kbb_a) * umask(i1:i2,j1:j2,jk) 
    555            END DO 
    556            tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hu_0(i1:i2,j1:j2) 
    557         END IF  
    558 # endif 
    559  
    560       ELSE 
    561  
    562 # if defined key_vertical 
    563          IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp 
    564  
    565          DO jj=j1,j2 
    566             DO ji=i1,i2 
    567                tabres_child(ji,jj,:) = 0._wp 
    568                N_in = mbku_parent(ji,jj) 
    569                zhtot = 0._wp 
    570                DO jk=1,N_in 
    571                   IF (jk==N_in) THEN 
    572                      h_in(jk) = hu0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot 
     657                  tabres_child(ji,jj,:) = 0._wp 
     658                  N_in = mbku_parent(ji,jj) 
     659                  zhtot = 0._wp 
     660                  DO jk=1,N_in 
     661                     !IF (jk==N_in) THEN 
     662                     !   h_in(jk) = hu0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot 
     663                     !ELSE 
     664                     !   h_in(jk) = tabres(ji,jj,jk,m2) 
     665                     !ENDIF 
     666                     h_in(jk) = e3u0_parent(ji,jj,jk) 
     667                     zhtot = zhtot + h_in(jk) 
     668                     tabin(jk) = tabres(ji,jj,jk,m1) 
     669                  END DO 
     670                  !          
     671                  N_out = 0 
     672                  DO jk=1,jpk 
     673                     IF (umask(ji,jj,jk) == 0) EXIT 
     674                     N_out = N_out + 1 
     675                     h_out(N_out) = e3u(ji,jj,jk,Kbb_a) 
     676                  END DO 
     677 
     678                  ! Account for small differences in free-surface 
     679                  IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 
     680                     h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 
    573681                  ELSE 
    574                      h_in(jk) = tabres(ji,jj,jk,m2) 
     682                     h_in(1)   = h_in(1)   - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 
    575683                  ENDIF 
    576                   zhtot = zhtot + h_in(jk) 
    577                   tabin(jk) = tabres(ji,jj,jk,m1) 
    578                END DO 
    579                !          
    580                N_out = 0 
    581                DO jk=1,jpk 
    582                   IF (umask(ji,jj,jk) == 0) EXIT 
    583                   N_out = N_out + 1 
    584                   h_out(N_out) = e3u(ji,jj,jk,Kbb_a) 
    585                END DO 
    586  
    587                ! Account for small differences in free-surface 
    588                IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 
    589                   h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 
    590                ELSE 
    591                   h_in(1)   = h_in(1)   - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 
    592                ENDIF 
    593684                   
    594                IF (N_in * N_out > 0) THEN 
    595                   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,1) 
    596                ENDIF  
    597             END DO 
    598          END DO 
    599  
    600          ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:) 
    601 #else 
    602          ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:) 
    603 #endif 
     685                  IF (N_in * N_out > 0) THEN 
     686                     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,1) 
     687                  ENDIF  
     688               END DO 
     689            END DO 
     690 
     691            ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:) 
     692         ELSE 
     693 
     694            ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:) 
     695   
     696         ENDIF 
    604697         ! 
    605698         DO jk = 1, jpkm1                                 ! Horizontal slab 
    … …  
    706799               DO ji=i1,i2 
    707800                  tabres(ji,jj,jk,m1) = vv(ji,jj,jk,Kbb_a) 
    708 # if defined key_vertical 
    709                   tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v(ji,jj,jk,Kbb_a) 
    710 # endif 
    711                END DO 
    712             END DO 
    713          END DO 
    714  
    715 # if defined key_vertical 
    716          ! Extrapolate thicknesses in partial bottom cells: 
    717          ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
    718          IF (ln_zps) THEN 
     801               END DO 
     802            END DO 
     803         END DO 
     804 
     805         IF ( l_vremap ) THEN 
     806 
     807            DO jk=k1,k2 
     808               DO jj=j1,j2 
     809                  DO ji=i1,i2 
     810                     tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v(ji,jj,jk,Kbb_a) 
     811                  END DO 
     812               END DO 
     813            END DO 
     814            ! Extrapolate thicknesses in partial bottom cells: 
     815            ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
     816            IF (ln_zps) THEN 
     817               DO jj=j1,j2 
     818                  DO ji=i1,i2 
     819                     jk = mbkv(ji,jj) 
     820                     tabres(ji,jj,jk,m2) = 0._wp 
     821                  END DO 
     822               END DO            
     823            END IF 
     824            ! Save ssh at last level: 
     825            tabres(i1:i2,j1:j2,k2,m2) = 0._wp 
     826            IF (.NOT.ln_linssh) THEN 
     827               ! This vertical sum below should be replaced by the sea-level at V-points (optimization): 
     828               DO jk=1,jpk 
     829                  tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3v(i1:i2,j1:j2,jk,Kbb_a) * vmask(i1:i2,j1:j2,jk) 
     830               END DO 
     831               tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hv_0(i1:i2,j1:j2) 
     832            END IF 
     833 
     834         END IF  
     835 
     836      ELSE 
     837 
     838         IF ( l_vremap ) THEN 
     839            IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp 
    719840            DO jj=j1,j2 
    720841               DO ji=i1,i2 
    721                   jk = mbkv(ji,jj) 
    722                   tabres(ji,jj,jk,m2) = 0._wp 
    723                END DO 
    724             END DO            
    725          END IF 
    726         ! Save ssh at last level: 
    727         tabres(i1:i2,j1:j2,k2,m2) = 0._wp 
    728         IF (.NOT.ln_linssh) THEN 
    729            ! This vertical sum below should be replaced by the sea-level at V-points (optimization): 
    730            DO jk=1,jpk 
    731               tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3v(i1:i2,j1:j2,jk,Kbb_a) * vmask(i1:i2,j1:j2,jk) 
    732            END DO 
    733            tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hv_0(i1:i2,j1:j2) 
    734         END IF  
    735 # endif 
    736  
    737       ELSE 
    738  
    739 # if defined key_vertical 
    740          IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp 
    741          DO jj=j1,j2 
    742             DO ji=i1,i2 
    743                tabres_child(ji,jj,:) = 0._wp 
    744                N_in = mbkv_parent(ji,jj) 
    745                zhtot = 0._wp 
    746                DO jk=1,N_in 
    747                   IF (jk==N_in) THEN 
    748                      h_in(jk) = hv0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot 
     842                  tabres_child(ji,jj,:) = 0._wp 
     843                  N_in = mbkv_parent(ji,jj) 
     844                  zhtot = 0._wp 
     845                  DO jk=1,N_in 
     846                     !IF (jk==N_in) THEN 
     847                     !   h_in(jk) = hv0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot 
     848                     !ELSE 
     849                     !   h_in(jk) = tabres(ji,jj,jk,m2) 
     850                     !ENDIF 
     851                     h_in(jk) = e3v0_parent(ji,jj,jk) 
     852                     zhtot = zhtot + h_in(jk) 
     853                     tabin(jk) = tabres(ji,jj,jk,m1) 
     854                  END DO 
     855                  !           
     856                  N_out = 0 
     857                  DO jk=1,jpk 
     858                     IF (vmask(ji,jj,jk) == 0) EXIT 
     859                     N_out = N_out + 1 
     860                     h_out(N_out) = e3v(ji,jj,jk,Kbb_a) 
     861                  END DO 
     862 
     863                  ! Account for small differences in free-surface 
     864                  IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 
     865                     h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 
    749866                  ELSE 
    750                      h_in(jk) = tabres(ji,jj,jk,m2) 
     867                     h_in(1)   = h_in(1) - (  sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 
    751868                  ENDIF 
    752                   zhtot = zhtot + h_in(jk) 
    753                   tabin(jk) = tabres(ji,jj,jk,m1) 
    754                END DO 
    755                !           
    756                N_out = 0 
    757                DO jk=1,jpk 
    758                   IF (vmask(ji,jj,jk) == 0) EXIT 
    759                   N_out = N_out + 1 
    760                   h_out(N_out) = e3v(ji,jj,jk,Kbb_a) 
    761                END DO 
    762  
    763                ! Account for small differences in free-surface 
    764                IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 
    765                   h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 
    766                ELSE 
    767                   h_in(1)   = h_in(1) - (  sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 
    768                ENDIF 
    769869          
    770                IF (N_in * N_out > 0) THEN 
    771                   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,1) 
    772                ENDIF 
    773             END DO 
    774          END DO 
    775  
    776          vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:)   
    777 # else 
    778          vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:) 
    779 # endif 
     870                  IF (N_in * N_out > 0) THEN 
     871                     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,1) 
     872                  ENDIF 
     873               END DO 
     874            END DO 
     875 
     876            vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:)   
     877         ELSE 
     878 
     879            vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:) 
     880 
     881         ENDIF 
    780882         ! 
    781883         DO jk = 1, jpkm1                                 ! Horizontal slab 
    … …  
    840942      ! 
    841943   END SUBROUTINE interpvn_sponge 
     944 
     945   SUBROUTINE interpunb_sponge(tabres,i1,i2,j1,j2, before) 
     946      !!--------------------------------------------- 
     947      !!   *** ROUTINE interpunb_sponge *** 
     948      !!---------------------------------------------     
     949      INTEGER, INTENT(in) :: i1,i2,j1,j2 
     950      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres 
     951      LOGICAL, INTENT(in) :: before 
     952 
     953      INTEGER  :: ji, jj, ind1, jmax 
     954      ! sponge parameters  
     955      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr 
     956      REAL(wp), DIMENSION(i1:i2,j1:j2) :: ubdiff 
     957      REAL(wp), DIMENSION(i1:i2,j1:j2) :: rotdiff, hdivdiff 
     958      !!---------------------------------------------     
     959      ! 
     960      IF( before ) THEN 
     961         DO jj=j1,j2 
     962            DO ji=i1,i2 
     963               tabres(ji,jj) = uu_b(ji,jj,Kmm_a) 
     964            END DO 
     965         END DO 
     966 
     967      ELSE 
     968 
     969         ubdiff(i1:i2,j1:j2) = (uu_b(i1:i2,j1:j2,Kmm_a) - tabres(i1:i2,j1:j2))*umask(i1:i2,j1:j2,1) 
     970         ! 
     971         !                                             ! -------- 
     972         ! Horizontal divergence                       !   div 
     973         !                                             ! -------- 
     974         DO jj = j1,j2 
     975            DO ji = i1+1,i2   ! vector opt. 
     976               zbtr = rn_sponge_dyn * r1_Dt * fspt_2d(ji,jj) * r1_ht_0(ji,jj) 
     977               hdivdiff(ji,jj) = (  e2u(ji  ,jj)*hu(ji  ,jj,Kbb_a) * ubdiff(ji  ,jj) & 
     978                                  &-e2u(ji-1,jj)*hu(ji-1,jj,Kbb_a) * ubdiff(ji-1,jj) ) * zbtr 
     979            END DO 
     980         END DO 
     981 
     982         DO jj = j1,j2-1 
     983            DO ji = i1,i2   ! vector opt. 
     984               zbtr = rn_sponge_dyn * r1_Dt * fspf_2d(ji,jj) * hf_0(ji,jj) 
     985               rotdiff(ji,jj) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1)   & 
     986                              &   +e1u(ji,jj  ) * ubdiff(ji,jj  ) ) * fmask(ji,jj,1) * zbtr  
     987            END DO 
     988         END DO 
     989         ! 
     990         DO jj = j1+1, j2-1 
     991            DO ji = i1+1, i2-1   ! vector opt. 
     992               IF (.NOT. tabspongedone_u(ji,jj)) THEN 
     993                  ze2u = rotdiff (ji,jj) 
     994                  ze1v = hdivdiff(ji,jj) 
     995                  ! horizontal diffusive trends 
     996                  zua = - ( ze2u - rotdiff (ji,jj-1) ) * r1_e2u(ji,jj) * r1_hu(ji,jj,Kmm_a)  & 
     997                      & + ( hdivdiff(ji+1,jj) - ze1v ) * r1_e1u(ji,jj)                       &  
     998                      & - rn_trelax_dyn * r1_Dt * fspu_2d(ji,jj) * ubdiff(ji,jj) 
     999 
     1000                  ! add it to the general momentum trends 
     1001                  uu(ji,jj,:,Krhs_a) = uu(ji,jj,:,Krhs_a) + zua                                  
     1002               ENDIF 
     1003            END DO 
     1004         END DO 
     1005 
     1006         tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE. 
     1007 
     1008         jmax = j2-1 
     1009         ind1 = jpjglo - ( nn_hls + nbghostcells + 2 )   ! North 
     1010         DO jj = mj0(ind1), mj1(ind1)                  
     1011            jmax = MIN(jmax,jj) 
     1012         END DO 
     1013 
     1014         DO jj = j1+1, jmax 
     1015            DO ji = i1+1, i2   ! vector opt. 
     1016               IF (.NOT. tabspongedone_v(ji,jj)) THEN 
     1017                     ze2u = rotdiff (ji,jj) 
     1018                     ze1v = hdivdiff(ji,jj) 
     1019                     zva = + ( ze2u - rotdiff (ji-1,jj) ) * r1_e1v(ji,jj) * r1_hv(ji,jj,Kmm_a) & 
     1020                           + ( hdivdiff(ji,jj+1) - ze1v ) * r1_e2v(ji,jj) 
     1021                     vv(ji,jj,:,Krhs_a) = vv(ji,jj,:,Krhs_a) + zva 
     1022               ENDIF 
     1023            END DO 
     1024         END DO 
     1025         ! 
     1026         tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE. 
     1027         ! 
     1028      ENDIF 
     1029      ! 
     1030   END SUBROUTINE interpunb_sponge 
     1031 
     1032    
     1033   SUBROUTINE interpvnb_sponge(tabres,i1,i2,j1,j2, before) 
     1034      !!--------------------------------------------- 
     1035      !!   *** ROUTINE interpvnb_sponge *** 
     1036      !!---------------------------------------------  
     1037      INTEGER, INTENT(in) :: i1,i2,j1,j2 
     1038      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres 
     1039      LOGICAL, INTENT(in) :: before 
     1040      ! 
     1041      INTEGER  ::   ji, jj, ind1, imax 
     1042      REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr 
     1043      REAL(wp), DIMENSION(i1:i2,j1:j2) :: vbdiff 
     1044      REAL(wp), DIMENSION(i1:i2,j1:j2) :: rotdiff, hdivdiff 
     1045      !!---------------------------------------------  
     1046       
     1047      IF( before ) THEN  
     1048         DO jj=j1,j2 
     1049            DO ji=i1,i2 
     1050               tabres(ji,jj) = vv_b(ji,jj,Kmm_a) 
     1051            END DO 
     1052         END DO 
     1053      ELSE 
     1054         vbdiff(i1:i2,j1:j2) = (vv_b(i1:i2,j1:j2,Kmm_a) - tabres(i1:i2,j1:j2))*vmask(i1:i2,j1:j2,1) 
     1055         !                                             ! -------- 
     1056         ! Horizontal divergence                       !   div 
     1057         !                                             ! -------- 
     1058         DO jj = j1+1,j2 
     1059            DO ji = i1,i2   ! vector opt. 
     1060               zbtr = rn_sponge_dyn * r1_Dt * fspt_2d(ji,jj) * r1_ht_0(ji,jj) 
     1061               hdivdiff(ji,jj) = ( e1v(ji,jj  ) * hv(ji,jj  ,Kbb_a) * vbdiff(ji,jj  )  & 
     1062                               &  -e1v(ji,jj-1) * hv(ji,jj-1,Kbb_a) * vbdiff(ji,jj-1)  ) * zbtr 
     1063            END DO 
     1064         END DO 
     1065         DO jj = j1,j2 
     1066            DO ji = i1,i2-1   ! vector opt. 
     1067               zbtr = rn_sponge_dyn * r1_Dt * fspf_2d(ji,jj) * hf_0(ji,jj)  
     1068               rotdiff(ji,jj) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj) &  
     1069                              &  -e2v(ji  ,jj) * vbdiff(ji  ,jj)  ) * fmask(ji,jj,1) * zbtr 
     1070            END DO 
     1071         END DO 
     1072         !                                                ! =============== 
     1073         !                                                 
     1074 
     1075         imax = i2 - 1 
     1076         ind1 = jpiglo - ( nn_hls + nbghostcells + 2 )   ! East 
     1077         DO ji = mi0(ind1), mi1(ind1)                 
     1078            imax = MIN(imax,ji) 
     1079         END DO 
     1080          
     1081         DO jj = j1+1, j2 
     1082            DO ji = i1+1, imax   ! vector opt. 
     1083               IF( .NOT. tabspongedone_u(ji,jj) ) THEN                                                      
     1084                  zua = - ( rotdiff (ji  ,jj) - rotdiff (ji,jj-1)) * r1_e2u(ji,jj) * r1_hu(ji,jj,Kmm_a)  & 
     1085                      & + ( hdivdiff(ji+1,jj) - hdivdiff(ji,jj  )) * r1_e1u(ji,jj) 
     1086                  uu(ji,jj,:,Krhs_a) = uu(ji,jj,:,Krhs_a) + zua 
     1087               ENDIF 
     1088            END DO 
     1089         END DO 
     1090         ! 
     1091         tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE. 
     1092         ! 
     1093         DO jj = j1+1, j2-1 
     1094            DO ji = i1+1, i2-1   ! vector opt. 
     1095               IF( .NOT. tabspongedone_v(ji,jj) ) THEN 
     1096                  zva  =  ( rotdiff (ji,jj  ) - rotdiff (ji-1,jj) ) * r1_e1v(ji,jj) *r1_hv(ji,jj,Kmm_a) & 
     1097                     &  + ( hdivdiff(ji,jj+1) - hdivdiff(ji  ,jj) ) * r1_e2v(ji,jj)                     & 
     1098                     &  - rn_trelax_dyn * r1_Dt * fspv_2d(ji,jj) * vbdiff(ji,jj) 
     1099                  vv(ji,jj,:,Krhs_a) = vv(ji,jj,:,Krhs_a) + zva 
     1100               ENDIF 
     1101            END DO 
     1102         END DO 
     1103         tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE. 
     1104      ENDIF 
     1105      ! 
     1106   END SUBROUTINE interpvnb_sponge 
     1107 
    8421108 
    8431109#else 

  • NEMO/trunk/src/NST/agrif_oce_update.F90

    r14064 r14086  
    1 #undef DECAL_FEEDBACK     /* SEPARATION of INTERFACES */ 
     1#undef DECAL_FEEDBACK    /* SEPARATION of INTERFACES */ 
    22#undef DECAL_FEEDBACK_2D  /* SEPARATION of INTERFACES (Barotropic mode) */ 
    33#undef VOL_REFLUX         /* VOLUME REFLUXING*/ 
    … …  
    2121   USE zdf_oce        ! vertical physics: ocean variables  
    2222   USE agrif_oce 
     23   USE dom_oce 
    2324   ! 
    2425   USE in_out_manager ! I/O manager 
    … …  
    3435 
    3536   PUBLIC   Agrif_Update_Tra, Agrif_Update_Dyn, Agrif_Update_vvl, Agrif_Update_ssh 
    36    PUBLIC   Update_Scales 
     37   PUBLIC   Update_Scales, Agrif_Check_parent_bat 
    3738 
    3839   !! * Substitutions 
    … …  
    5455      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update tracers  from grid Number',Agrif_Fixed() 
    5556 
    56 #if defined key_vertical 
    57 ! Effect of this has to be carrefully checked  
    58 ! depending on what the nesting tools ensure for 
    59 ! volume conservation: 
     57      l_vremap                      = ln_vert_remap 
     58      Agrif_UseSpecialValueInUpdate = .NOT.l_vremap 
     59      Agrif_SpecialValueFineGrid    = 0._wp 
     60      !  
     61# if ! defined DECAL_FEEDBACK 
     62      CALL Agrif_Update_Variable(ts_update_id, procname=updateTS) 
     63! near boundary update: 
     64!      CALL Agrif_Update_Variable(ts_update_id,locupdate=(/0,2/), procname=updateTS) 
     65# else 
     66      CALL Agrif_Update_Variable(ts_update_id, locupdate=(/1,0/),procname=updateTS) 
     67! near boundary update: 
     68!      CALL Agrif_Update_Variable(ts_update_id,locupdate=(/1,2/), procname=updateTS) 
     69# endif 
     70      ! 
    6071      Agrif_UseSpecialValueInUpdate = .FALSE. 
    61 #else 
    62       Agrif_UseSpecialValueInUpdate = .TRUE. 
    63 #endif 
    64       Agrif_SpecialValueFineGrid    = 0._wp 
    65       !  
    66 # if ! defined DECAL_FEEDBACK 
    67       CALL Agrif_Update_Variable(tsn_id, procname=updateTS) 
    68 ! near boundary update: 
    69 !      CALL Agrif_Update_Variable(tsn_id,locupdate=(/0,2/), procname=updateTS) 
    70 # else 
    71       CALL Agrif_Update_Variable(tsn_id, locupdate=(/1,0/),procname=updateTS) 
    72 ! near boundary update: 
    73 !      CALL Agrif_Update_Variable(tsn_id,locupdate=(/1,2/), procname=updateTS) 
    74 # endif 
    75       ! 
    76       Agrif_UseSpecialValueInUpdate = .FALSE. 
     72      l_vremap                      = .FALSE. 
    7773      ! 
    7874      ! 
    … …  
    9086      Agrif_UseSpecialValueInUpdate = .FALSE. 
    9187      Agrif_SpecialValueFineGrid    = 0._wp 
    92  
    93       use_sign_north = .TRUE. 
    94       sign_north     = -1._wp 
    95  
    96       !      
     88      l_vremap                      = ln_vert_remap 
     89      use_sign_north                = .TRUE. 
     90      sign_north                    = -1._wp      
     91! 
     92# if ! defined DECAL_FEEDBACK_2D 
     93      CALL Agrif_Update_Variable(unb_update_id,locupdate1=(/  nn_shift_bar,-2/),locupdate2=(/  nn_shift_bar,-2/),procname = updateU2d) 
     94      CALL Agrif_Update_Variable(vnb_update_id,locupdate1=(/  nn_shift_bar,-2/),locupdate2=(/  nn_shift_bar,-2/),procname = updateV2d) 
     95# else 
     96      CALL Agrif_Update_Variable(unb_update_id,locupdate1=(/  nn_shift_bar,-2/),locupdate2=(/1+nn_shift_bar,-2/),procname = updateU2d) 
     97      CALL Agrif_Update_Variable(vnb_update_id,locupdate1=(/1+nn_shift_bar,-2/),locupdate2=(/  nn_shift_bar,-2/),procname = updateV2d)   
     98# endif 
     99      !  
     100      IF ( ln_dynspg_ts .AND. ln_bt_fw ) THEN 
     101         ! Update time integrated transports 
     102#  if ! defined DECAL_FEEDBACK_2D 
     103         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/  nn_shift_bar,-2/),locupdate2=(/  nn_shift_bar,-2/),procname = updateub2b) 
     104         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/  nn_shift_bar,-2/),locupdate2=(/  nn_shift_bar,-2/),procname = updatevb2b) 
     105#  else 
     106         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/  nn_shift_bar,-2/),locupdate2=(/1+nn_shift_bar,-2/),procname = updateub2b) 
     107         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1+nn_shift_bar,-2/),locupdate2=(/  nn_shift_bar,-2/),procname = updatevb2b) 
     108#  endif 
     109      END IF 
     110 
    97111# if ! defined DECAL_FEEDBACK 
    98112      CALL Agrif_Update_Variable(un_update_id,procname = updateU) 
    … …  
    108122!      CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updateV) 
    109123# endif 
    110  
    111 # if ! defined DECAL_FEEDBACK_2D 
    112       CALL Agrif_Update_Variable(e1u_id,procname = updateU2d) 
    113       CALL Agrif_Update_Variable(e2v_id,procname = updateV2d)   
    114 # else 
    115       CALL Agrif_Update_Variable(e1u_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU2d) 
    116       CALL Agrif_Update_Variable(e2v_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV2d)   
    117 # endif 
    118       ! 
    119 # if ! defined DECAL_FEEDBACK_2D 
    120       ! Account for updated thicknesses at boundary edges 
    121       IF (.NOT.ln_linssh) THEN 
    122 !         CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_u_bdy) 
    123 !         CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_v_bdy) 
    124       ENDIF 
    125 # endif 
    126       !  
    127       IF ( ln_dynspg_ts .AND. ln_bt_fw ) THEN 
    128          ! Update time integrated transports 
    129 #  if ! defined DECAL_FEEDBACK_2D 
    130          CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b) 
    131          CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b) 
    132 #  else 
    133          CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateub2b) 
    134          CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updatevb2b) 
    135 #  endif 
    136       END IF 
    137124      ! 
    138125      use_sign_north = .FALSE. 
     126      l_vremap = .FALSE. 
    139127      ! 
    140128   END SUBROUTINE Agrif_Update_Dyn 
    … …  
    150138      Agrif_SpecialValueFineGrid = 0._wp 
    151139# if ! defined DECAL_FEEDBACK_2D 
    152       CALL Agrif_Update_Variable(sshn_id,procname = updateSSH) 
     140      CALL Agrif_Update_Variable(sshn_id,locupdate=(/  nn_shift_bar,-2/), procname = updateSSH)  
    153141# else 
    154       CALL Agrif_Update_Variable(sshn_id,locupdate=(/1,0/),procname = updateSSH) 
     142      CALL Agrif_Update_Variable(sshn_id,locupdate=(/1+nn_shift_bar,-2/),procname = updateSSH) 
    155143# endif 
    156144      ! 
    … …  
    158146      ! 
    159147#  if defined VOL_REFLUX 
    160       IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN 
     148      IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN  
    161149         use_sign_north = .TRUE. 
    162150         sign_north = -1._wp 
    163151         ! Refluxing on ssh: 
    164152#  if defined DECAL_FEEDBACK_2D 
    165          CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0, 0/),locupdate2=(/1, 1/),procname = reflux_sshu) 
    166          CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1, 1/),locupdate2=(/0, 0/),procname = reflux_sshv) 
     153         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/nn_shift_bar,nn_shift_bar/),locupdate2=(/1+nn_shift_bar,1+nn_shift_bar/),procname = reflux_sshu) 
     154         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1+nn_shift_bar,1+nn_shift_bar/),locupdate2=(/nn_shift_bar,nn_shift_bar/),procname = reflux_sshv) 
    167155#  else 
    168          CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/-1,-1/),locupdate2=(/ 0, 0/),procname = reflux_sshu) 
    169          CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/ 0, 0/),locupdate2=(/-1,-1/),procname = reflux_sshv) 
     156         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/-1+nn_shift_bar,-1+nn_shift_bar/),locupdate2=(/nn_shift_bar, nn_shift_bar/),procname = reflux_sshu) 
     157         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/ nn_shift_bar, nn_shift_bar/),locupdate2=(/-1+nn_shift_bar,-1+nn_shift_bar/),procname = reflux_sshv) 
    170158#  endif 
    171159         use_sign_north = .FALSE. 
    … …  
    320308#endif 
    321309 
    322 #if defined key_vertical 
    323310 
    324311   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    … …  
    340327      ! 
    341328      IF (before) THEN 
    342 !jc_alt 
    343 !         AGRIF_SpecialValue = -999._wp 
    344          DO jn = n1,n2-1 
     329         IF ( l_vremap ) THEN 
     330            DO jn = n1,n2-1 
     331               DO jk=k1,k2 
     332                  DO jj=j1,j2 
     333                     DO ji=i1,i2 
     334                        tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) 
     335                     END DO 
     336                  END DO 
     337               END DO 
     338            END DO 
    345339            DO jk=k1,k2 
    346340               DO jj=j1,j2 
    347341                  DO ji=i1,i2 
    348 !jc_alt 
    349 !                     tabres(ji,jj,jk,jn) = (ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) ) & 
    350 !                                         &  * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1._wp) * 999._wp 
    351                      tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) 
     342                     tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 
    352343                  END DO 
    353344               END DO 
    354345            END DO 
    355          END DO 
    356          DO jk=k1,k2 
     346         ELSE 
     347            DO jn = 1,jpts 
     348               DO jk=k1,k2 
     349                  DO jj=j1,j2 
     350                     DO ji=i1,i2 
     351                        tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a)  * e3t(ji,jj,jk,Kmm_a) / e3t_0(ji,jj,jk) 
     352                     END DO 
     353                  END DO 
     354               END DO 
     355            END DO 
     356 
     357         ENDIF 
     358      ELSE 
     359         IF ( l_vremap ) THEN 
     360            tabres_child(:,:,:,:) = 0._wp 
     361            AGRIF_SpecialValue = 0._wp 
    357362            DO jj=j1,j2 
    358363               DO ji=i1,i2 
    359 !jc_alt 
    360 !                  tabres(ji,jj,jk,n2) =      tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) & 
    361 !                                      &   + (tmask(ji,jj,jk) - 1._wp) * 999._wp 
    362                   tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 
    363                END DO 
    364             END DO 
    365          END DO 
    366       ELSE 
    367          tabres_child(:,:,:,:) = 0._wp 
    368          AGRIF_SpecialValue = 0._wp 
    369          DO jj=j1,j2 
    370             DO ji=i1,i2 
    371                N_in = 0 
    372                DO jk=k1,k2 !k2 = jpk of child grid 
    373 ! jc_alt 
    374 !                  IF (tabres(ji,jj,jk,n2) < -900._wp  ) EXIT 
    375                   IF (tabres(ji,jj,jk,n2) == 0._wp  ) EXIT 
    376                   N_in = N_in + 1 
    377                   tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 
    378                   h_in(N_in) = tabres(ji,jj,jk,n2) 
     364                  N_in = 0 
     365                  DO jk=k1,k2 !k2 = jpk of child grid 
     366                     IF (tabres(ji,jj,jk,n2) <= 1.e-6_wp  ) EXIT 
     367                     N_in = N_in + 1 
     368                     tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 
     369                     h_in(N_in) = tabres(ji,jj,jk,n2) 
     370                  ENDDO 
     371                  N_out = 0 
     372                  DO jk=1,jpk ! jpk of parent grid 
     373                     IF (tmask(ji,jj,jk) == 0 ) EXIT ! TODO: Will not work with ISF 
     374                     N_out = N_out + 1 
     375                     h_out(N_out) = e3t(ji,jj,jk,Kmm_a)  
     376                  ENDDO 
     377                  IF (N_in*N_out > 0) THEN !Remove this? 
     378                     CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts) 
     379                  ENDIF 
    379380               ENDDO 
    380                N_out = 0 
    381                DO jk=1,jpk ! jpk of parent grid 
    382                   IF (tmask(ji,jj,jk) == 0 ) EXIT ! TODO: Will not work with ISF 
    383                   N_out = N_out + 1 
    384                   h_out(N_out) = e3t(ji,jj,jk,Kmm_a)  
    385                ENDDO 
    386                IF (N_in*N_out > 0) THEN !Remove this? 
    387                   CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts) 
    388                ENDIF 
    389381            ENDDO 
    390          ENDDO 
    391  
    392          IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN 
    393             ! Add asselin part 
     382 
     383            IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN 
     384               ! Add asselin part 
     385               DO jn = 1,jpts 
     386                  DO jk = 1, jpkm1 
     387                     DO jj = j1, j2 
     388                        DO ji = i1, i2 
     389                           IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN 
     390                              ztb  = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 
     391                              ztnu = tabres_child(ji,jj,jk,jn) * e3t(ji,jj,jk,Kmm_a) 
     392                              ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a) 
     393                              ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) )  &  
     394                                        &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a) 
     395                           ENDIF 
     396                        END DO 
     397                     END DO 
     398                  END DO 
     399               END DO 
     400            ENDIF 
    394401            DO jn = 1,jpts 
    395402               DO jk = 1, jpkm1 
    396403                  DO jj = j1, j2 
    397404                     DO ji = i1, i2 
    398                         IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN 
    399                            ztb  = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 
    400                            ztnu = tabres_child(ji,jj,jk,jn) * e3t(ji,jj,jk,Kmm_a) 
    401                            ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a) 
    402                            ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) )  &  
    403                                      &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a) 
    404                         ENDIF 
     405                        IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN  
     406                           ts(ji,jj,jk,jn,Kmm_a) = tabres_child(ji,jj,jk,jn) 
     407                        END IF 
    405408                     END DO 
    406409                  END DO 
    407410               END DO 
    408411            END DO 
    409          ENDIF 
    410          DO jn = 1,jpts 
    411             DO jk = 1, jpkm1 
    412                DO jj = j1, j2 
    413                   DO ji = i1, i2 
    414                      IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN  
    415                         ts(ji,jj,jk,jn,Kmm_a) = tabres_child(ji,jj,jk,jn) 
    416                      END IF 
     412         ELSE 
     413            DO jn = 1,jpts 
     414               tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) & 
     415                                            & * tmask(i1:i2,j1:j2,k1:k2) 
     416            ENDDO 
     417  
     418            IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN 
     419               ! Add asselin part 
     420               DO jn = 1,jpts 
     421                  DO jk = k1, k2 
     422                     DO jj = j1, j2 
     423                        DO ji = i1, i2 
     424                           IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN 
     425                              ztb  = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 
     426                              ztnu = tabres(ji,jj,jk,jn) 
     427                              ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a) 
     428                              ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) )  &  
     429                                        &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a) 
     430                           ENDIF 
     431                        END DO 
     432                     END DO 
    417433                  END DO 
    418434               END DO 
    419             END DO 
    420          END DO 
    421          ! 
     435            ENDIF 
     436            DO jn = 1,jpts 
     437               DO jk=k1,k2 
     438                  DO jj=j1,j2 
     439                     DO ji=i1,i2 
     440                        IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN  
     441                           ts(ji,jj,jk,jn,Kmm_a) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm_a) 
     442                        END IF 
     443                     END DO 
     444                  END DO 
     445               END DO 
     446            END DO 
     447            ! 
     448         ENDIF 
    422449         IF  ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 
    423450            ts(i1:i2,j1:j2,1:jpkm1,1:jpts,Kbb_a)  = ts(i1:i2,j1:j2,1:jpkm1,1:jpts,Kmm_a) 
    424          ENDIF 
     451         ENDIF          
    425452      ENDIF 
    426453      !  
    427454   END SUBROUTINE updateTS 
    428455 
    429 # else 
    430  
    431    SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    432       !!--------------------------------------------- 
    433       !!           *** ROUTINE updateT *** 
    434       !!--------------------------------------------- 
    435       INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2 
    436       REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    437       LOGICAL, INTENT(in) :: before 
    438       !! 
    439       INTEGER :: ji,jj,jk,jn 
    440       REAL(wp) :: ztb, ztnu, ztno 
    441       !!--------------------------------------------- 
    442       ! 
    443       IF (before) THEN 
    444          DO jn = 1,jpts 
    445             DO jk=k1,k2 
    446                DO jj=j1,j2 
    447                   DO ji=i1,i2 
    448 !> jc tmp 
    449                      tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a)  * e3t(ji,jj,jk,Kmm_a) / e3t_0(ji,jj,jk) 
    450 !                     tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a)  * e3t(ji,jj,jk,Kmm_a) 
    451 !< jc tmp 
    452                   END DO 
    453                END DO 
    454             END DO 
    455          END DO 
    456       ELSE 
    457 !> jc tmp 
    458          DO jn = 1,jpts 
    459             tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) & 
    460                                          & * tmask(i1:i2,j1:j2,k1:k2) 
    461          ENDDO 
    462 !< jc tmp 
    463          IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN 
    464             ! Add asselin part 
    465             DO jn = 1,jpts 
    466                DO jk = k1, k2 
    467                   DO jj = j1, j2 
    468                      DO ji = i1, i2 
    469                         IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN 
    470                            ztb  = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 
    471                            ztnu = tabres(ji,jj,jk,jn) 
    472                            ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a) 
    473                            ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) )  &  
    474                                      &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a) 
    475                         ENDIF 
    476                      END DO 
    477                   END DO 
    478                END DO 
    479             END DO 
    480          ENDIF 
    481          DO jn = 1,jpts 
    482             DO jk=k1,k2 
    483                DO jj=j1,j2 
    484                   DO ji=i1,i2 
    485                      IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN  
    486                         ts(ji,jj,jk,jn,Kmm_a) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm_a) 
    487                      END IF 
    488                   END DO 
    489                END DO 
    490             END DO 
    491          END DO 
    492          ! 
    493          IF  ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 
    494             ts(i1:i2,j1:j2,k1:k2,1:jpts,Kbb_a)  = ts(i1:i2,j1:j2,k1:k2,1:jpts,Kmm_a) 
    495          ENDIF 
    496          ! 
    497       ENDIF 
    498       !  
    499    END SUBROUTINE updateTS 
    500  
    501 # endif 
    502  
    503 # if defined key_vertical 
    504456 
    505457   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    … …  
    513465      INTEGER ::   ji, jj, jk 
    514466      REAL(wp)::   zrhoy, zub, zunu, zuno 
     467      REAL(wp), DIMENSION(jpi,jpj) ::   zpgu   ! 2D workspace 
    515468! VERTICAL REFINEMENT BEGIN 
    516469      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 
    517470      REAL(wp) :: h_in(k1:k2) 
    518471      REAL(wp) :: h_out(1:jpk) 
    519       INTEGER  :: N_in, N_out 
    520       REAL(wp) :: h_diff, excess, thick 
     472      INTEGER  :: N_in, N_out, N_in_save, N_out_save 
     473      REAL(wp) :: zhmin, zd 
    521474      REAL(wp) :: tabin(k1:k2) 
    522475! VERTICAL REFINEMENT END 
    … …  
    525478      IF( before ) THEN 
    526479         zrhoy = Agrif_Rhoy() 
    527 !jc_alt 
    528 !         AGRIF_SpecialValue = -999._wp 
    529480         DO jk=k1,k2 
     481            tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm_a) &  
     482                                     &   * umask(i1:i2,j1:j2,jk) * uu(i1:i2,j1:j2,jk,Kmm_a)   
     483         END DO 
     484 
     485         IF ( l_vremap ) THEN 
     486            DO jk=k1,k2 
     487               tabres(i1:i2,j1:j2,jk,2) = zrhoy * umask(i1:i2,j1:j2,jk) * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm_a) 
     488            END DO 
     489         ENDIF 
     490 
     491      ELSE 
     492 
     493         tabres_child(:,:,:) = 0._wp 
     494 
     495         IF ( l_vremap ) THEN 
     496 
    530497            DO jj=j1,j2 
    531498               DO ji=i1,i2 
    532 !jc_alt 
    533 !                  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)  & 
    534 !                                     &  + (umask(ji,jj,jk)-1._wp)*999._wp 
    535                   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)   
    536 !jc_alt 
    537 !                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)  & 
    538 !                                     &  + (umask(ji,jj,jk)-1._wp)*999._wp 
    539                   tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) 
    540                END DO 
    541             END DO 
    542          END DO 
    543       ELSE 
    544          tabres_child(:,:,:) = 0. 
    545          AGRIF_SpecialValue = 0._wp 
    546          DO jj=j1,j2 
    547             DO ji=i1,i2 
    548                N_in = 0 
    549                h_in(:) = 0._wp 
    550                tabin(:) = 0._wp 
    551                DO jk=k1,k2 !k2=jpk of child grid 
    552 !jc_alt 
    553 !                  IF( tabres(ji,jj,jk,2) < -900._wp) EXIT 
    554                   IF( tabres(ji,jj,jk,2) == 0.) EXIT 
    555                   N_in = N_in + 1 
    556                   tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
    557                   h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj) 
     499                  N_in = 0 
     500                  h_in(:) = 0._wp 
     501                  tabin(:) = 0._wp 
     502                  DO jk=k1,k2 !k2=jpk of child grid 
     503                     IF( tabres(ji,jj,jk,2)*r1_e2u(ji,jj) <= 1.e-6_wp ) EXIT 
     504                     N_in = N_in + 1 
     505                     tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
     506                     h_in(N_in) = tabres(ji,jj,jk,2) * r1_e2u(ji,jj) 
     507                  ENDDO 
     508                  N_out = 0 
     509                  DO jk=1,jpk 
     510                     IF (umask(ji,jj,jk) == 0._wp) EXIT 
     511                     N_out = N_out + 1 
     512                     h_out(N_out) = e3u(ji,jj,jk,Kmm_a) 
     513                  ENDDO 
     514                  IF (N_in * N_out > 0) THEN 
     515                     ! Deal with potentially different depths at velocity points: 
     516                     N_in_save  = N_in 
     517                     N_out_save = N_out 
     518                     IF ( ABS(sum(h_out(1:N_out))-sum(h_in(1:N_in))) > 1.e-6_wp ) THEN 
     519                        zhmin = MIN(sum(h_out(1:N_out)), sum(h_in(1:N_in))) 
     520                        zd = 0._wp 
     521                        DO jk=1, N_in_save 
     522                           IF ( (zd +  h_in(jk)) > zhmin-1.e-6) THEN 
     523                              N_in = jk 
     524                              h_in(jk) = zhmin - zd 
     525                              EXIT  
     526                           ENDIF 
     527                           zd = zd + h_in(jk) 
     528                        END DO 
     529                        zd = 0._wp 
     530                        DO jk=1, N_out_save 
     531                           IF ( (zd +  h_out(jk)) > zhmin-1.e-6) THEN 
     532                              N_out = jk 
     533                              h_out(jk) = zhmin - zd 
     534                              EXIT  
     535                           ENDIF 
     536                           zd = zd + h_out(jk) 
     537                        END DO 
     538                     END IF 
     539                     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,1) 
     540                     IF (N_out < N_out_save) tabres_child(ji,jj,N_out+1:N_out_save) = tabres_child(ji,jj,N_out) 
     541                  ENDIF 
    558542               ENDDO 
    559                N_out = 0 
    560                DO jk=1,jpk 
    561                   IF (umask(ji,jj,jk) == 0) EXIT 
    562                   N_out = N_out + 1 
    563                   h_out(N_out) = e3u(ji,jj,jk,Kmm_a) 
    564                ENDDO 
    565                IF (N_in * N_out > 0) THEN 
    566                   h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    567                   excess = 0._wp 
    568                   IF (h_diff < -1.e-4) THEN 
    569 !Even if bathy at T points match it's possible for the U points to be deeper in the child grid.  
    570 !In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. 
    571                      DO jk=N_in,1,-1 
    572                         thick = MIN(-1*h_diff, h_in(jk)) 
    573                         excess = excess + tabin(jk)*thick*e2u(ji,jj) 
    574                         tabin(jk) = tabin(jk)*(1. - thick/h_in(jk)) 
    575                         h_diff = h_diff + thick 
    576                         IF ( h_diff == 0) THEN 
    577                            N_in = jk 
    578                            h_in(jk) = h_in(jk) - thick 
    579                            EXIT 
    580                         ENDIF 
    581                      ENDDO 
    582                   ENDIF 
    583                   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,1) 
    584                   tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out)) 
    585                ENDIF 
    586543            ENDDO 
    587          ENDDO 
     544 
     545         ELSE 
     546            DO jk=1,jpk 
     547               DO jj=j1,j2 
     548                  DO ji=i1,i2 
     549                     tabres_child(ji,jj,jk) = tabres(ji,jj,jk,1) * r1_e2u(ji,jj) /  e3u(ji,jj,jk,Kmm_a) 
     550                  END DO 
     551               END DO 
     552            END DO 
     553         ENDIF 
    588554         ! 
    589555         DO jk=1,jpk 
    … …  
    603569         END DO 
    604570         ! 
     571         ! Correct now and before transports: 
     572         DO jj=j1,j2 
     573            DO ji=i1,i2 
     574               zpgu(ji,jj) = 0._wp 
     575               DO jk=1,jpkm1 
     576                  zpgu(ji,jj) = zpgu(ji,jj) + e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a) 
     577               END DO 
     578               ! 
     579               DO jk=1,jpkm1               
     580                  uu(ji,jj,jk,Kmm_a) = uu(ji,jj,jk,Kmm_a) + & 
     581                      &  (uu_b(ji,jj,Kmm_a) - zpgu(ji,jj) * r1_hu(ji,jj,Kmm_a)) * umask(ji,jj,jk)            
     582               END DO 
     583               ! 
     584               zpgu(ji,jj) = 0._wp 
     585               DO jk=1,jpkm1 
     586                  zpgu(ji,jj) = zpgu(ji,jj) + e3u(ji,jj,jk,Kbb_a) * uu(ji,jj,jk,Kbb_a) 
     587               END DO 
     588               ! 
     589               DO jk=1,jpkm1               
     590                  uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) + & 
     591                      &  (uu_b(ji,jj,Kbb_a) - zpgu(ji,jj) * r1_hu(ji,jj,Kbb_a)) * umask(ji,jj,jk)            
     592               END DO 
     593               ! 
     594            END DO 
     595         END DO 
     596         ! 
    605597         IF  ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 
    606598            uu(i1:i2,j1:j2,1:jpkm1,Kbb_a)  = uu(i1:i2,j1:j2,1:jpkm1,Kmm_a) 
    … …  
    611603   END SUBROUTINE updateu 
    612604 
    613 #else 
    614  
    615    SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    616       !!--------------------------------------------- 
    617       !!           *** ROUTINE updateu *** 
    618       !!--------------------------------------------- 
    619       INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
    620       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    621       LOGICAL                                     , INTENT(in   ) :: before 
    622       ! 
    623       INTEGER  :: ji, jj, jk 
    624       REAL(wp) :: zrhoy, zub, zunu, zuno 
    625       !!--------------------------------------------- 
    626       !  
    627       IF( before ) THEN 
    628          zrhoy = Agrif_Rhoy() 
    629          DO jk = k1, k2 
    630             tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm_a) * uu(i1:i2,j1:j2,jk,Kmm_a) 
    631          END DO 
    632       ELSE 
    633          DO jk=k1,k2 
    634             DO jj=j1,j2 
    635                DO ji=i1,i2 
    636                   tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e2u(ji,jj)  
    637                   ! 
    638                   IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part 
    639                      zub  = uu(ji,jj,jk,Kbb_a) * e3u(ji,jj,jk,Kbb_a)  ! fse3t_b prior update should be used 
    640                      zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a) 
    641                      zunu = tabres(ji,jj,jk,1) 
    642                      uu(ji,jj,jk,Kbb_a) = ( zub + rn_atfp * ( zunu - zuno) ) &       
    643                                     & * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a) 
    644                   ENDIF 
    645                   ! 
    646                   uu(ji,jj,jk,Kmm_a) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u(ji,jj,jk,Kmm_a) 
    647                END DO 
    648             END DO 
    649          END DO 
    650          ! 
    651          IF  ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 
    652             uu(i1:i2,j1:j2,k1:k2,Kbb_a)  = uu(i1:i2,j1:j2,k1:k2,Kmm_a) 
    653          ENDIF 
    654          ! 
    655       ENDIF 
    656       !  
    657    END SUBROUTINE updateu 
    658  
    659 # endif 
    660  
    661    SUBROUTINE correct_u_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) 
    662       !!--------------------------------------------- 
    663       !!           *** ROUTINE correct_u_bdy *** 
     605 
     606   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
     607      !!--------------------------------------------- 
     608      !!           *** ROUTINE updatev *** 
    664609      !!--------------------------------------------- 
    665610      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
    666611      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    667612      LOGICAL                                     , INTENT(in   ) :: before 
    668       INTEGER                                     , INTENT(in)    :: nb, ndir 
    669       !! 
    670       LOGICAL :: western_side, eastern_side  
    671       ! 
    672       INTEGER  :: jj, jk 
    673       REAL(wp) :: zcor 
    674       !!--------------------------------------------- 
    675       !  
    676       IF( .NOT.before ) THEN 
    677          ! 
    678          western_side  = (nb == 1).AND.(ndir == 1) 
    679          eastern_side  = (nb == 1).AND.(ndir == 2) 
    680          ! 
    681          IF (western_side) THEN 
    682             DO jj=j1,j2 
    683                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) 
    684                uu_b(i1-1,jj,Kmm_a) = uu_b(i1-1,jj,Kmm_a) + zcor 
    685                DO jk=1,jpkm1 
    686                   uu(i1-1,jj,jk,Kmm_a) = uu(i1-1,jj,jk,Kmm_a) + zcor * umask(i1-1,jj,jk) 
    687                END DO  
    688             END DO 
    689          ENDIF 
    690          ! 
    691          IF (eastern_side) THEN 
    692             DO jj=j1,j2 
    693                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) 
    694                uu_b(i2+1,jj,Kmm_a) = uu_b(i2+1,jj,Kmm_a) + zcor 
    695                DO jk=1,jpkm1 
    696                   uu(i2+1,jj,jk,Kmm_a) = uu(i2+1,jj,jk,Kmm_a) + zcor * umask(i2+1,jj,jk) 
    697                END DO  
    698             END DO 
    699          ENDIF 
    700          ! 
    701       ENDIF 
    702       !  
    703    END SUBROUTINE correct_u_bdy 
    704  
    705 # if  defined key_vertical 
    706  
    707    SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    708       !!--------------------------------------------- 
    709       !!           *** ROUTINE updatev *** 
    710       !!--------------------------------------------- 
    711       INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
    712       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    713       LOGICAL                                     , INTENT(in   ) :: before 
    714613      ! 
    715614      INTEGER  ::   ji, jj, jk 
    716615      REAL(wp) ::   zrhox, zvb, zvnu, zvno 
     616      REAL(wp), DIMENSION(jpi,jpj) ::   zpgv   ! 2D workspace 
    717617! VERTICAL REFINEMENT BEGIN 
    718618      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 
    719619      REAL(wp) :: h_in(k1:k2) 
    720620      REAL(wp) :: h_out(1:jpk) 
    721       INTEGER :: N_in, N_out 
    722       REAL(wp) :: h_diff, excess, thick 
     621      INTEGER  :: N_in, N_out, N_in_save, N_out_save 
     622      REAL(wp) :: zhmin, zd 
    723623      REAL(wp) :: tabin(k1:k2) 
    724624! VERTICAL REFINEMENT END 
    … …  
    727627      IF( before ) THEN 
    728628         zrhox = Agrif_Rhox() 
    729 !jc_alt 
    730 !         AGRIF_SpecialValue = -999._wp 
    731629         DO jk=k1,k2 
     630            tabres(i1:i2,j1:j2,jk,1) = zrhox * e1v(i1:i2,j1:j2) * e3v(i1:i2,j1:j2,jk,Kmm_a) &  
     631                                     &   * vmask(i1:i2,j1:j2,jk) * vv(i1:i2,j1:j2,jk,Kmm_a)   
     632         END DO 
     633 
     634         IF ( l_vremap ) THEN 
     635            DO jk=k1,k2 
     636               tabres(i1:i2,j1:j2,jk,2) = zrhox * e1v(i1:i2,j1:j2) * e3v(i1:i2,j1:j2,jk,Kmm_a) * vmask(i1:i2,j1:j2,jk) 
     637            END DO 
     638         ENDIF 
     639 
     640      ELSE 
     641 
     642         tabres_child(:,:,:) = 0._wp 
     643 
     644         IF ( l_vremap ) THEN 
     645 
    732646            DO jj=j1,j2 
    733647               DO ji=i1,i2 
    734 !jc_alt 
    735 !                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) * vv(ji,jj,jk,Kmm_a) & 
    736 !                                     & + (vmask(ji,jj,jk)-1._wp) * 999._wp 
    737                   tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) * vv(ji,jj,jk,Kmm_a)  
    738 !jc_alt 
    739 !                  tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) & 
    740 !                                     & + (vmask(ji,jj,jk)-1._wp) * 999._wp 
    741                   tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) 
    742                END DO 
    743             END DO 
    744          END DO 
    745       ELSE 
    746          tabres_child(:,:,:) = 0. 
    747          AGRIF_SpecialValue = 0._wp 
    748          DO jj=j1,j2 
    749             DO ji=i1,i2 
    750                N_in = 0 
    751                DO jk=k1,k2 
    752 !jc_alt 
    753 !                  IF (tabres(ji,jj,jk,2) < -900._wp) EXIT 
    754                   IF (tabres(ji,jj,jk,2) == 0) EXIT 
    755                   N_in = N_in + 1 
    756                   tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
    757                   h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj) 
     648                  N_in = 0 
     649                  DO jk=k1,k2 
     650                     IF (tabres(ji,jj,jk,2)* r1_e1v(ji,jj) <= 1.e-6_wp) EXIT 
     651                     N_in = N_in + 1 
     652                     tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
     653                     h_in(N_in) = tabres(ji,jj,jk,2) * r1_e1v(ji,jj) 
     654                  ENDDO 
     655                  N_out = 0 
     656                  DO jk=1,jpk 
     657                     IF (vmask(ji,jj,jk) == 0) EXIT 
     658                     N_out = N_out + 1 
     659                     h_out(N_out) = e3v(ji,jj,jk,Kmm_a) 
     660                  ENDDO 
     661                  IF (N_in * N_out > 0) THEN 
     662                     ! Deal with potentially different depths at velocity points: 
     663                     N_in_save  = N_in 
     664                     N_out_save = N_out 
     665                     IF ( ABS(sum(h_out(1:N_out))-sum(h_in(1:N_in))) > 1.e-6_wp ) THEN 
     666                        zhmin = MIN(sum(h_out(1:N_out)), sum(h_in(1:N_in))) 
     667                        zd = 0._wp 
     668                        DO jk=1, N_in_save 
     669                           IF ( (zd +  h_in(jk)) > zhmin-1.e-6) THEN 
     670                              N_in = jk 
     671                              h_in(jk) = zhmin - zd 
     672                              EXIT  
     673                           ENDIF 
     674                           zd = zd + h_in(jk) 
     675                        END DO 
     676                        zd = 0._wp 
     677                        DO jk=1, N_out_save 
     678                           IF ( (zd +  h_out(jk)) > zhmin-1.e-6) THEN 
     679                              N_out = jk 
     680                              h_out(jk) = zhmin - zd 
     681                              EXIT  
     682                           ENDIF 
     683                           zd = zd + h_out(jk) 
     684                        END DO 
     685                     END IF 
     686                     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,1) 
     687                     IF (N_out < N_out_save) tabres_child(ji,jj,N_out+1:N_out_save) = tabres_child(ji,jj,N_out) 
     688                  ENDIF 
    758689               ENDDO 
    759                N_out = 0 
    760                DO jk=1,jpk 
    761                   IF (vmask(ji,jj,jk) == 0) EXIT 
    762                   N_out = N_out + 1 
    763                   h_out(N_out) = e3v(ji,jj,jk,Kmm_a) 
    764                ENDDO 
    765                IF (N_in * N_out > 0) THEN 
    766                   h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    767                   excess = 0._wp 
    768                   IF (h_diff < -1.e-4) then 
    769 !Even if bathy at T points match it's possible for the V points to be deeper in the child grid.  
    770 !In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. 
    771                      DO jk=N_in,1,-1 
    772                         thick = MIN(-1*h_diff, h_in(jk)) 
    773                         excess = excess + tabin(jk)*thick*e2u(ji,jj) 
    774                         tabin(jk) = tabin(jk)*(1. - thick/h_in(jk)) 
    775                         h_diff = h_diff + thick 
    776                         IF ( h_diff == 0) THEN 
    777                            N_in = jk 
    778                            h_in(jk) = h_in(jk) - thick 
    779                            EXIT 
    780                         ENDIF 
    781                      ENDDO 
    782                   ENDIF 
    783                   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,1) 
    784                   tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out)) 
    785                ENDIF 
    786690            ENDDO 
    787          ENDDO 
     691 
     692         ELSE 
     693            DO jk=1,jpk 
     694               DO jj=j1,j2 
     695                  DO ji=i1,i2 
     696                     tabres_child(ji,jj,jk) = tabres(ji,jj,jk,1) * r1_e1v(ji,jj) /  e3v(ji,jj,jk,Kmm_a) 
     697                  END DO 
     698               END DO 
     699            END DO 
     700         ENDIF 
    788701         ! 
    789702         DO jk=1,jpkm1 
    … …  
    803716         END DO 
    804717         ! 
     718         ! Correct now and before transports: 
     719         DO jj=j1,j2 
     720            DO ji=i1,i2 
     721               zpgv(ji,jj) = 0._wp 
     722               DO jk=1,jpkm1 
     723                  zpgv(ji,jj) = zpgv(ji,jj) + e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a) 
     724               END DO 
     725               ! 
     726               DO jk=1,jpkm1               
     727                  vv(ji,jj,jk,Kmm_a) = vv(ji,jj,jk,Kmm_a) + & 
     728                      &  (vv_b(ji,jj,Kmm_a) - zpgv(ji,jj) * r1_hv(ji,jj,Kmm_a)) * vmask(ji,jj,jk)            
     729               END DO 
     730               ! 
     731               zpgv(ji,jj) = 0._wp 
     732               DO jk=1,jpkm1 
     733                  zpgv(ji,jj) = zpgv(ji,jj) + e3v(ji,jj,jk,Kbb_a) * vv(ji,jj,jk,Kbb_a) 
     734               END DO 
     735               ! 
     736               DO jk=1,jpkm1               
     737                  vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) + & 
     738                      &  (vv_b(ji,jj,Kbb_a) - zpgv(ji,jj) * r1_hv(ji,jj,Kbb_a)) * vmask(ji,jj,jk)            
     739               END DO 
     740               ! 
     741            END DO 
     742         END DO 
     743         ! 
    805744         IF  ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 
    806745            vv(i1:i2,j1:j2,1:jpkm1,Kbb_a)  = vv(i1:i2,j1:j2,1:jpkm1,Kmm_a) 
    … …  
    810749      !  
    811750   END SUBROUTINE updatev 
    812  
    813 # else 
    814  
    815    SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before) 
    816       !!--------------------------------------------- 
    817       !!           *** ROUTINE updatev *** 
    818       !!--------------------------------------------- 
    819       INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
    820       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    821       LOGICAL                                     , INTENT(in   ) :: before 
    822       ! 
    823       INTEGER  :: ji, jj, jk 
    824       REAL(wp) :: zrhox, zvb, zvnu, zvno 
    825       !!---------------------------------------------       
    826       ! 
    827       IF (before) THEN 
    828          zrhox = Agrif_Rhox() 
    829          DO jk=k1,k2 
    830             DO jj=j1,j2 
    831                DO ji=i1,i2 
    832                   tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a) 
    833                END DO 
    834             END DO 
    835          END DO 
    836       ELSE 
    837          DO jk=k1,k2 
    838             DO jj=j1,j2 
    839                DO ji=i1,i2 
    840                   tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e1v(ji,jj) 
    841                   ! 
    842                   IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part 
    843                      zvb  = vv(ji,jj,jk,Kbb_a) * e3v(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 
    844                      zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a) 
    845                      zvnu = tabres(ji,jj,jk,1) 
    846                      vv(ji,jj,jk,Kbb_a) = ( zvb + rn_atfp * ( zvnu - zvno) ) &       
    847                                     & * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a) 
    848                   ENDIF 
    849                   ! 
    850                   vv(ji,jj,jk,Kmm_a) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kmm_a) 
    851                END DO 
    852             END DO 
    853          END DO 
    854          ! 
    855          IF  ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 
    856             vv(i1:i2,j1:j2,k1:k2,Kbb_a)  = vv(i1:i2,j1:j2,k1:k2,Kmm_a) 
    857          ENDIF 
    858          ! 
    859       ENDIF 
    860       !  
    861    END SUBROUTINE updatev 
    862  
    863 # endif 
    864  
    865    SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) 
    866       !!--------------------------------------------- 
    867       !!           *** ROUTINE correct_v_bdy *** 
    868       !!--------------------------------------------- 
    869       INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
    870       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    871       LOGICAL                                     , INTENT(in   ) :: before 
    872       INTEGER                                     , INTENT(in)    :: nb, ndir 
    873       !! 
    874       LOGICAL :: southern_side, northern_side  
    875       ! 
    876       INTEGER  :: ji, jk 
    877       REAL(wp) :: zcor 
    878       !!--------------------------------------------- 
    879       !  
    880       IF( .NOT.before ) THEN 
    881          ! 
    882          southern_side = (nb == 2).AND.(ndir == 1) 
    883          northern_side = (nb == 2).AND.(ndir == 2) 
    884          ! 
    885          IF (southern_side) THEN 
    886             DO ji=i1,i2 
    887                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) 
    888                vv_b(ji,j1-1,Kmm_a) = vv_b(ji,j1-1,Kmm_a) + zcor 
    889                DO jk=1,jpkm1 
    890                   vv(ji,j1-1,jk,Kmm_a) = vv(ji,j1-1,jk,Kmm_a) + zcor * vmask(ji,j1-1,jk) 
    891                END DO  
    892             END DO 
    893          ENDIF 
    894          ! 
    895          IF (northern_side) THEN 
    896             DO ji=i1,i2 
    897                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) 
    898                vv_b(ji,j2+1,Kmm_a) = vv_b(ji,j2+1,Kmm_a) + zcor 
    899                DO jk=1,jpkm1 
    900                   vv(ji,j2+1,jk,Kmm_a) = vv(ji,j2+1,jk,Kmm_a) + zcor * vmask(ji,j2+1,jk) 
    901                END DO  
    902             END DO 
    903          ENDIF 
    904          ! 
    905       ENDIF 
    906       !  
    907    END SUBROUTINE correct_v_bdy 
    908751 
    909752 
    … …  
    935778               tabres(ji,jj) =  tabres(ji,jj) * r1_e2u(ji,jj)   
    936779               !     
    937                ! Update "now" 3d velocities: 
    938                zpgu(ji,jj) = 0._wp 
    939                DO jk=1,jpkm1 
    940                   zpgu(ji,jj) = zpgu(ji,jj) + e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a) 
    941                END DO 
    942                ! 
    943                zcorr = (tabres(ji,jj) - zpgu(ji,jj)) * r1_hu(ji,jj,Kmm_a) 
    944                DO jk=1,jpkm1               
    945                   uu(ji,jj,jk,Kmm_a) = uu(ji,jj,jk,Kmm_a) + zcorr * umask(ji,jj,jk)            
    946                END DO 
    947                ! 
    948780               ! Update barotropic velocities: 
    949781               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 
    … …  
    955787               uu_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hu(ji,jj,Kmm_a) * umask(ji,jj,1) 
    956788               !        
    957                ! Correct "before" velocities to hold correct bt component: 
    958                zpgu(ji,jj) = 0.e0 
    959                DO jk=1,jpkm1 
    960                   zpgu(ji,jj) = zpgu(ji,jj) + e3u(ji,jj,jk,Kbb_a) * uu(ji,jj,jk,Kbb_a) 
    961                END DO 
    962                ! 
    963                zcorr = uu_b(ji,jj,Kbb_a) - zpgu(ji,jj) * r1_hu(ji,jj,Kbb_a) 
    964                DO jk=1,jpkm1               
    965                   uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) + zcorr * umask(ji,jj,jk)            
    966                END DO 
    967                ! 
    968789            END DO 
    969790         END DO 
    … …  
    1002823            DO ji=i1,i2 
    1003824               tabres(ji,jj) =  tabres(ji,jj) * r1_e1v(ji,jj)   
    1004                !     
    1005                ! Update "now" 3d velocities: 
    1006                zpgv(ji,jj) = 0.e0 
    1007                DO jk=1,jpkm1 
    1008                   zpgv(ji,jj) = zpgv(ji,jj) + e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a) 
    1009                END DO 
    1010                ! 
    1011                zcorr = (tabres(ji,jj) - zpgv(ji,jj)) * r1_hv(ji,jj,Kmm_a) 
    1012                DO jk=1,jpkm1               
    1013                   vv(ji,jj,jk,Kmm_a) = vv(ji,jj,jk,Kmm_a) + zcorr * vmask(ji,jj,jk)            
    1014                END DO 
    1015                ! 
    1016825               ! Update barotropic velocities: 
    1017826               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 
    … …  
    1023832               vv_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hv(ji,jj,Kmm_a) * vmask(ji,jj,1) 
    1024833               !        
    1025                ! Correct "before" velocities to hold correct bt component: 
    1026                zpgv(ji,jj) = 0.e0 
    1027                DO jk=1,jpkm1 
    1028                   zpgv(ji,jj) = zpgv(ji,jj) + e3v(ji,jj,jk,Kbb_a) * vv(ji,jj,jk,Kbb_a) 
    1029                END DO 
    1030                ! 
    1031                zcorr = vv_b(ji,jj,Kbb_a) - zpgv(ji,jj) * r1_hv(ji,jj,Kbb_a) 
    1032                DO jk=1,jpkm1               
    1033                   vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) + zcorr * vmask(ji,jj,jk)            
    1034                END DO 
    1035                ! 
    1036834            END DO 
    1037835         END DO 
    … …  
    1120918               ub2_i_b(ji,jj) = ub2_i_b(ji,jj) + za1 * zcor  
    1121919               ! Update corrective fluxes: 
    1122                un_bf(ji,jj)  = un_bf(ji,jj) + zcor 
     920               IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) un_bf(ji,jj)  = un_bf(ji,jj) + zcor 
    1123921               ! Update half step back fluxes: 
    1124922               ub2_b(ji,jj) = tabres(ji,jj) 
    … …  
    12081006               vb2_i_b(ji,jj) = vb2_i_b(ji,jj) + za1 * zcor  
    12091007               ! Update corrective fluxes: 
    1210                vn_bf(ji,jj)  = vn_bf(ji,jj) + zcor 
     1008               IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler)))  vn_bf(ji,jj)  = vn_bf(ji,jj) + zcor 
    12111009               ! Update half step back fluxes: 
    12121010               vb2_b(ji,jj) = tabres(ji,jj) 
    … …  
    14791277#endif 
    14801278 
     1279   SUBROUTINE Agrif_Check_parent_bat( ) 
     1280      !!---------------------------------------------------------------------- 
     1281      !!                   *** ROUTINE Agrif_Check_parent_bat *** 
     1282      !!---------------------------------------------------------------------- 
     1283      !  
     1284      IF (( .NOT.ln_agrif_2way ).OR.(.NOT.ln_chk_bathy).OR.(Agrif_Root())) RETURN 
     1285      ! 
     1286      Agrif_UseSpecialValueInUpdate = .FALSE. 
     1287      ! 
     1288      IF(lwp) WRITE(numout,*) ' ' 
     1289      IF(lwp) WRITE(numout,*) 'AGRIF: Check parent volume at Level:', Agrif_Level() 
     1290      ! 
     1291# if ! defined DECAL_FEEDBACK 
     1292      CALL Agrif_Update_Variable(batupd_id,procname = update_bat) 
     1293# else 
     1294      CALL Agrif_Update_Variable(batupd_id,locupdate=(/1,0/),procname = update_bat) 
     1295# endif 
     1296      ! 
     1297      kindic_agr = Agrif_Parent(kindic_agr) 
     1298      CALL mpp_sum( 'Agrif_Check_parent_bat', kindic_agr ) 
     1299 
     1300      IF( kindic_agr /= 0 ) THEN 
     1301         CALL ctl_stop('==> Averaged Bathymetry does not match parent volume')  
     1302      ELSE 
     1303         IF(lwp) WRITE(numout,*) '==> Averaged Bathymetry matches parent '  
     1304         IF(lwp) WRITE(numout,*) '' 
     1305      ENDIF 
     1306      ! 
     1307   END SUBROUTINE Agrif_Check_parent_bat 
     1308 
     1309   SUBROUTINE update_bat(ptab, i1, i2, j1, j2, before ) 
     1310      !!--------------------------------------------- 
     1311      !!           *** ROUTINE update_bat *** 
     1312      !!--------------------------------------------- 
     1313      REAL(wp), DIMENSION(i1:i2,j1:j2) :: ptab 
     1314      INTEGER, INTENT(in) :: i1, i2, j1, j2 
     1315      LOGICAL, INTENT(in) :: before 
     1316      INTEGER :: ji, jj 
     1317      ! 
     1318      !!--------------------------------------------- 
     1319      ! 
     1320      IF( before ) THEN 
     1321         ptab(i1:i2,j1:j2) = ht_0(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1) 
     1322      ELSE 
     1323         kindic_agr = 0 
     1324         ! 
     1325         DO jj=j1,j2 
     1326            DO ji=i1,i2 
     1327               IF ( (ssmask(ji,jj).NE.0._wp).AND.& 
     1328               & (ABS(ptab(ji,jj)-ht_0(ji,jj)).GE.1.e-6) ) THEN  
     1329                  kindic_agr = kindic_agr + 1  
     1330               ENDIF 
     1331            END DO 
     1332         END DO 
     1333         ! 
     1334      ENDIF 
     1335      ! 
     1336   END SUBROUTINE update_bat 
     1337 
    14811338#else 
    14821339   !!---------------------------------------------------------------------- 

  • NEMO/trunk/src/NST/agrif_top_interp.F90

    r13216 r14086  
    4343      Agrif_SpecialValue    = 0._wp 
    4444      Agrif_UseSpecialValue = .TRUE. 
     45      l_vremap              = ln_vert_remap 
    4546      ! 
    4647      CALL Agrif_Bc_variable( trn_id, procname=interptrn ) 
     48      ! 
    4749      Agrif_UseSpecialValue = .FALSE. 
     50      l_vremap              = .FALSE. 
    4851      ! 
    4952   END SUBROUTINE Agrif_trc 
    … …  
    5760      LOGICAL                                     , INTENT(in   ) ::   before 
    5861      ! 
    59       INTEGER  ::   ji, jj, jk, jn, ibdy, jbdy   ! dummy loop indices 
    60       INTEGER  ::   imin, imax, jmin, jmax, N_in, N_out 
    61       REAL(wp) ::   zrho, z1, z2, z3, z4, z5, z6, z7 
    62  
     62      INTEGER  ::   ji, jj, jk, jn  ! dummy loop indices 
     63      INTEGER  ::   N_in, N_out 
     64      INTEGER  :: item 
    6365      ! vertical interpolation: 
    64       REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,1:jptra) :: ptab_child 
    65       REAL(wp), DIMENSION(k1:k2,1:jptra) :: tabin 
    66       REAL(wp), DIMENSION(k1:k2) :: h_in 
    67       REAL(wp), DIMENSION(1:jpk) :: h_out 
    68       !!---------------------------------------------------------------------- 
    69  
    70       IF( before ) THEN          
     66      REAL(wp) :: zhtot, zwgt 
     67      REAL(wp), DIMENSION(k1:k2,1:jptra) :: tabin, tabin_i 
     68      REAL(wp), DIMENSION(k1:k2) :: z_in, h_in_i, z_in_i 
     69      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out 
     70      !!---------------------------------------------------------------------- 
     71 
     72      IF( before ) THEN 
     73 
     74         item = Kmm_a 
     75         IF( l_ini_child )   Kmm_a = Kbb_a   
     76 
    7177         DO jn = 1,jptra 
    7278            DO jk=k1,k2 
    … …  
    7783              END DO 
    7884           END DO 
    79         END DO 
    80  
    81 # if defined key_vertical 
    82         DO jk=k1,k2 
    83            DO jj=j1,j2 
    84               DO ji=i1,i2 
    85                  ptab(ji,jj,jk,jptra+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a)  
    86               END DO 
    87            END DO 
    88         END DO 
    89 # endif 
     85         END DO 
     86 
     87         IF( l_vremap .OR. l_ini_child .OR. ln_zps ) THEN 
     88            ! Fill cell depths (i.e. gdept) to be interpolated 
     89            ! Warning: these are masked, hence extrapolated prior interpolation. 
     90            DO jj=j1,j2 
     91               DO ji=i1,i2 
     92                  ptab(ji,jj,k1,jptra+1) = 0.5_wp * tmask(ji,jj,k1) * e3t(ji,jj,k1,Kmm_a) 
     93                  DO jk=k1+1,k2 
     94                     ptab(ji,jj,jk,jptra+1) = tmask(ji,jj,jk) * & 
     95                        & ( ptab(ji,jj,jk-1,jptra+1) + 0.5_wp * (e3t(ji,jj,jk-1,Kmm_a)+e3t(ji,jj,jk,Kmm_a)) ) 
     96                  END DO 
     97               END DO 
     98            END DO 
     99 
     100            ! Save ssh at last level: 
     101            IF (.NOT.ln_linssh) THEN 
     102               ptab(i1:i2,j1:j2,k2,jptra+1) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1)  
     103            END IF       
     104         ENDIF 
     105         Kmm_a = item 
     106 
    90107      ELSE  
    91  
    92 # if defined key_vertical 
    93          DO jj=j1,j2 
    94             DO ji=i1,i2 
    95                ptab_child(ji,jj,:) = 0._wp 
    96                N_in = 0 
    97                DO jk=k1,k2 !k2 = jpk of parent grid 
    98                   IF (ptab(ji,jj,jk,n2) == 0) EXIT 
    99                   N_in = N_in + 1 
    100                   tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1) 
    101                   h_in(N_in) = ptab(ji,jj,jk,n2) 
     108         item = Krhs_a 
     109         IF( l_ini_child )   Krhs_a = Kbb_a   
     110 
     111         IF( l_vremap .OR. l_ini_child ) THEN 
     112            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,n2) = 0._wp  
     113                
     114            DO jj=j1,j2 
     115               DO ji=i1,i2 
     116                  tr(ji,jj,:,:,Krhs_a) = 0.   
     117                  ! 
     118                  ! Build vertical grids: 
     119                  N_in = mbkt_parent(ji,jj) 
     120                  ! Input grid (account for partial cells if any): 
     121                  DO jk=1,N_in 
     122                     z_in(jk) = ptab(ji,jj,jk,n2) - ptab(ji,jj,k2,n2) 
     123                     tabin(jk,1:jptra) = ptab(ji,jj,jk,1:jptra) 
     124                  END DO 
     125                   
     126                  ! Intermediate grid: 
     127                  IF ( l_vremap ) THEN 
     128                     DO jk = 1, N_in 
     129                        h_in_i(jk) = e3t0_parent(ji,jj,jk) * &  
     130                          &       (1._wp + ptab(ji,jj,k2,n2)/(ht0_parent(ji,jj)*ssmask(ji,jj) + 1._wp - ssmask(ji,jj))) 
     131                     END DO 
     132                     z_in_i(1) = 0.5_wp * h_in_i(1) 
     133                     DO jk=2,N_in 
     134                        z_in_i(jk) = z_in_i(jk-1) + 0.5_wp * ( h_in_i(jk) + h_in_i(jk-1) ) 
     135                     END DO 
     136                     z_in_i(1:N_in) = z_in_i(1:N_in)  - ptab(ji,jj,k2,n2) 
     137                  ENDIF 
     138 
     139                  ! Output (Child) grid: 
     140                  N_out = mbkt(ji,jj) 
     141                  DO jk=1,N_out 
     142                     h_out(jk) = e3t(ji,jj,jk,Krhs_a) 
     143                  END DO 
     144                  z_out(1) = 0.5_wp * h_out(1) 
     145                  DO jk=2,N_out 
     146                     z_out(jk) = z_out(jk-1) + 0.5_wp * ( h_out(jk)+h_out(jk-1) ) 
     147                  END DO 
     148                  IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out)  - ssh(ji,jj,Krhs_a)                
     149 
     150                  IF (N_in*N_out > 0) THEN 
     151                     IF( l_ini_child ) THEN 
     152                        CALL remap_linear(tabin(1:N_in,1:jptra),z_in(1:N_in),tr(ji,jj,1:N_out,1:jptra,Krhs_a),          & 
     153                                      &   z_out(1:N_out),N_in,N_out,jptra)   
     154                     ELSE   
     155                        CALL remap_linear(tabin(1:N_in,1:jptra),z_in(1:N_in),tabin_i(1:N_in,1:jptra),                     & 
     156                                     &   z_in_i(1:N_in),N_in,N_in,jptra) 
     157                        CALL reconstructandremap(tabin_i(1:N_in,1:jptra),h_in_i(1:N_in),tr(ji,jj,1:N_out,1:jptra,Krhs_a), & 
     158                                      &   h_out(1:N_out),N_in,N_out,jptra)    
     159                     ENDIF 
     160                  ENDIF 
    102161               END DO 
    103                N_out = 0 
    104                DO jk=1,jpk ! jpk of child grid 
    105                   IF (tmask(ji,jj,jk) == 0) EXIT  
    106                   N_out = N_out + 1 
    107                   h_out(jk) = e3t(ji,jj,jk,Krhs_a) 
    108                ENDDO 
    109                IF (N_in > 0) THEN 
    110                   CALL reconstructandremap(tabin(1:N_in,1:jptra),h_in,ptab_child(ji,jj,1:N_out,1:jptra),h_out,N_in,N_out,jptra) 
    111                ENDIF 
    112             ENDDO 
    113          ENDDO 
    114 # else 
    115          ptab_child(i1:i2,j1:j2,1:jpk,1:jptra) = ptab(i1:i2,j1:j2,1:jpk,1:jptra) 
    116 # endif 
    117          ! 
    118          DO jn=1, jptra 
    119             tr(i1:i2,j1:j2,1:jpk,jn,Krhs_a)=ptab_child(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk)  
    120          END DO 
     162            END DO 
     163            Krhs_a = item 
     164  
     165         ELSE 
     166 
     167            IF ( Agrif_Parent(ln_zps) ) THEN ! Account for partial cells  
     168                                             ! linear vertical interpolation 
     169               DO jj=j1,j2 
     170                  DO ji=i1,i2 
     171                     ! 
     172                     N_in  = mbkt(ji,jj) 
     173                     N_out = mbkt(ji,jj) 
     174                     z_in(1) = ptab(ji,jj,1,n2) 
     175                     tabin(1,1:jptra) = ptab(ji,jj,1,1:jptra) 
     176                     DO jk=2, N_in 
     177                        z_in(jk) = ptab(ji,jj,jk,n2) 
     178                        tabin(jk,1:jptra) = ptab(ji,jj,jk,1:jptra) 
     179                     END DO 
     180                     IF (.NOT.ln_linssh) z_in(1:N_in) = z_in(1:N_in) - ptab(ji,jj,k2,n2) 
     181                     z_out(1) = 0.5_wp * e3t(ji,jj,1,Krhs_a) 
     182                     DO jk=2, N_out 
     183                        z_out(jk) = z_out(jk-1) + 0.5_wp * (e3t(ji,jj,jk-1,Krhs_a) + e3t(ji,jj,jk,Krhs_a)) 
     184                     END DO 
     185                     IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out) - ssh(ji,jj,Krhs_a) 
     186                     CALL remap_linear(tabin(1:N_in,1:jptra),z_in(1:N_in),ptab(ji,jj,1:N_out,1:jptra), & 
     187                                   &   z_out(1:N_out),N_in,N_out,jptra)   
     188                  END DO 
     189               END DO 
     190 
     191            ENDIF 
     192 
     193            DO jn=1, jptra 
     194                tr(i1:i2,j1:j2,1:jpk,jn,Krhs_a)=ptab(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk)  
     195            END DO 
     196         ENDIF 
     197 
    121198      ENDIF 
    122199      ! 

  • NEMO/trunk/src/NST/agrif_top_sponge.F90

    r12489 r14086  
    4545      ! 
    4646#if defined SPONGE_TOP 
    47 !! Assume persistence  
     47      !! Assume persistence: 
    4848      zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot()) 
    49       CALL Agrif_sponge 
     49 
    5050      Agrif_SpecialValue    = 0._wp 
    5151      Agrif_UseSpecialValue = .TRUE. 
     52      l_vremap              = ln_vert_remap 
    5253      tabspongedone_trn     = .FALSE. 
     54      ! 
    5355      CALL Agrif_Bc_Variable( trn_sponge_id, calledweight=zcoef, procname=interptrn_sponge ) 
     56      ! 
    5457      Agrif_UseSpecialValue = .FALSE. 
     58      l_vremap              = .FALSE. 
    5559#endif 
    5660      ! 
    … …  
    5862 
    5963 
    60    SUBROUTINE interptrn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    61       !!---------------------------------------------------------------------- 
    62       !!                   *** ROUTINE interptrn_sponge *** 
     64   SUBROUTINE interptrn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before)  
     65      !!---------------------------------------------------------------------- 
     66      !!                 *** ROUTINE interptrn_sponge *** 
    6367      !!---------------------------------------------------------------------- 
    6468      INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2 
    … …  
    6771      ! 
    6872      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    69       REAL(wp) ::   zabe1, zabe2, ztrelax 
    70       REAL(wp), DIMENSION(i1:i2,j1:j2)               ::   ztu, ztv 
    71       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,1:jptra) ::   trbdiff 
     73      INTEGER  ::   iku, ikv 
     74      REAL(wp) :: ztra, zabe1, zabe2, zbtr, zhtot 
     75      REAL(wp), DIMENSION(i1-1:i2,j1-1:j2,jpk) :: ztu, ztv 
     76      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::trbdiff 
    7277      ! vertical interpolation: 
    73       REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,1:jptra) ::tabres_child 
    74       REAL(wp), DIMENSION(k1:k2,1:jptra) :: tabin 
    75       REAL(wp), DIMENSION(k1:k2) :: h_in 
    76       REAL(wp), DIMENSION(1:jpk) :: h_out 
     78      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child 
     79      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin, tabin_i 
     80      REAL(wp), DIMENSION(k1:k2) :: z_in, z_in_i, h_in_i 
     81      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out 
    7782      INTEGER :: N_in, N_out 
    78       REAL(wp) :: h_diff 
    7983      !!---------------------------------------------------------------------- 
    8084      ! 
    … …  
    9094         END DO 
    9195 
    92 # if defined key_vertical 
    93          DO jk=k1,k2 
    94             DO jj=j1,j2 
    95                DO ji=i1,i2 
    96                   tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kbb_a)  
    97                END DO 
    98             END DO 
    99          END DO 
    100 # endif 
    101       ELSE       
    102 # if defined key_vertical 
    103          tabres_child(:,:,:,:) = 0. 
    104          DO jj=j1,j2 
    105             DO ji=i1,i2 
    106                N_in = 0 
    107                DO jk=k1,k2 !k2 = jpk of parent grid 
    108                   IF (tabres(ji,jj,jk,n2) == 0) EXIT 
    109                   N_in = N_in + 1 
    110                   tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1) 
    111                   h_in(N_in) = tabres(ji,jj,jk,n2) 
    112                END DO 
    113                N_out = 0 
    114                DO jk=1,jpk ! jpk of child grid 
    115                   IF (tmask(ji,jj,jk) == 0) EXIT  
    116                   N_out = N_out + 1 
    117                   h_out(jk) = e3t(ji,jj,jk,Kbb_a) !Child grid scale factors. Could multiply by e1e2t here instead of division above 
    118                ENDDO 
    119                IF (N_in > 0) THEN 
    120                   CALL reconstructandremap(tabin(1:N_in,1:jptra),h_in,tabres_child(ji,jj,1:N_out,1:jptra),h_out,N_in,N_out,jptra) 
     96         IF ( l_vremap.OR.ln_zps ) THEN 
     97 
     98            ! Fill cell depths (i.e. gdept) to be interpolated 
     99            ! Warning: these are masked, hence extrapolated prior interpolation. 
     100            DO jj=j1,j2 
     101               DO ji=i1,i2 
     102                  tabres(ji,jj,k1,jptra+1) = 0.5_wp * tmask(ji,jj,k1) * e3t(ji,jj,k1,Kbb_a) 
     103                  DO jk=k1+1,k2 
     104                     tabres(ji,jj,jk,jptra+1) = tmask(ji,jj,jk) * & 
     105                        & ( tabres(ji,jj,jk-1,jptra+1) + 0.5_wp * (e3t(ji,jj,jk-1,Kbb_a)+e3t(ji,jj,jk,Kbb_a)) ) 
     106                  END DO 
     107               END DO 
     108            END DO 
     109 
     110            ! Save ssh at last level: 
     111            IF ( .NOT.ln_linssh ) THEN 
     112               tabres(i1:i2,j1:j2,k2,jptra+1) = ssh(i1:i2,j1:j2,Kbb_a)*tmask(i1:i2,j1:j2,1)  
     113            END IF   
     114     
     115         END IF 
     116 
     117      ELSE  
     118         ! 
     119         IF ( l_vremap ) THEN 
     120 
     121            IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp 
     122 
     123            DO jj=j1,j2 
     124               DO ji=i1,i2 
     125 
     126                  tabres_child(ji,jj,:,:) = 0._wp  
     127                  ! Build vertical grids: 
     128                  N_in = mbkt_parent(ji,jj) 
     129                  ! Input grid (account for partial cells if any): 
     130                  DO jk=1,N_in 
     131                     z_in(jk) = tabres(ji,jj,jk,n2) - tabres(ji,jj,k2,n2) 
     132                     tabin(jk,1:jptra) = tabres(ji,jj,jk,1:jptra) 
     133                  END DO 
     134                   
     135                  ! Intermediate grid: 
     136                  DO jk = 1, N_in 
     137                     h_in_i(jk) = e3t0_parent(ji,jj,jk) * &  
     138                       &       (1._wp + tabres(ji,jj,k2,n2)/(ht0_parent(ji,jj)*ssmask(ji,jj) + 1._wp - ssmask(ji,jj))) 
     139                  END DO 
     140                  z_in_i(1) = 0.5_wp * h_in_i(1) 
     141                  DO jk=2,N_in 
     142                     z_in_i(jk) = z_in_i(jk-1) + 0.5_wp * ( h_in_i(jk) + h_in_i(jk-1) ) 
     143                  END DO 
     144                  z_in_i(1:N_in) = z_in_i(1:N_in)  - tabres(ji,jj,k2,n2) 
     145 
     146                  ! Output (Child) grid: 
     147                  N_out = mbkt(ji,jj) 
     148                  DO jk=1,N_out 
     149                     h_out(jk) = e3t(ji,jj,jk,Kbb_a) 
     150                  END DO 
     151                  z_out(1) = 0.5_wp * h_out(1) 
     152                  DO jk=2,N_out 
     153                     z_out(jk) = z_out(jk-1) + 0.5_wp * ( h_out(jk)+h_out(jk-1) ) 
     154                  END DO 
     155                  IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out)  - ssh(ji,jj,Kbb_a) 
     156 
     157                  ! Account for small differences in the free-surface 
     158                  IF ( sum(h_out(1:N_out)) > sum(h_in_i(1:N_in) )) THEN 
     159                     h_out(1) = h_out(1)  - ( sum(h_out(1:N_out))-sum(h_in_i(1:N_in)) ) 
     160                  ELSE 
     161                     h_in_i(1)= h_in_i(1) - ( sum(h_in_i(1:N_in))-sum(h_out(1:N_out)) ) 
     162                  END IF 
     163                  IF (N_in*N_out > 0) THEN 
     164                     CALL remap_linear(tabin(1:N_in,1:jptra),z_in(1:N_in),tabin_i(1:N_in,1:jptra),z_in_i(1:N_in),N_in,N_in,jptra) 
     165                     CALL reconstructandremap(tabin_i(1:N_in,1:jptra),h_in_i(1:N_in),tabres_child(ji,jj,1:N_out,1:jptra),h_out(1:N_out),N_in,N_out,jptra) 
     166!                     CALL remap_linear(tabin(1:N_in,1:jptra),z_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jptra),z_out(1:N_in),N_in,N_out,jptra)   
     167                  ENDIF 
     168               END DO 
     169            END DO 
     170 
     171            DO jj=j1,j2 
     172               DO ji=i1,i2 
     173                  DO jk=1,jpkm1 
     174                     trbdiff(ji,jj,jk,1:jptra) = (tr(ji,jj,jk,1:jptra,Kbb_a) - tabres_child(ji,jj,jk,1:jptra)) * tmask(ji,jj,jk) 
     175                  END DO 
     176               END DO 
     177            END DO 
     178 
     179         ELSE 
     180 
     181            IF ( Agrif_Parent(ln_zps) ) THEN ! Account for partial cells 
     182 
     183               DO jj=j1,j2 
     184                  DO ji=i1,i2 
     185                     ! 
     186                     N_in  = mbkt(ji,jj)  
     187                     N_out = mbkt(ji,jj)  
     188                     z_in(1) = tabres(ji,jj,1,n2) 
     189                     tabin(1,1:jptra) = tabres(ji,jj,1,1:jptra) 
     190                     DO jk=2, N_in 
     191                        z_in(jk) = tabres(ji,jj,jk,n2) 
     192                        tabin(jk,1:jptra) = tabres(ji,jj,jk,1:jptra) 
     193                     END DO  
     194                     IF (.NOT.ln_linssh) z_in(1:N_in) = z_in(1:N_in) - tabres(ji,jj,k2,n2) 
     195 
     196                     z_out(1) = 0.5_wp * e3t(ji,jj,1,Kbb_a) 
     197                     DO jk=2, N_out 
     198                        z_out(jk) = z_out(jk-1) + 0.5_wp * (e3t(ji,jj,jk-1,Kbb_a) + e3t(ji,jj,jk,Kbb_a))  
     199                     END DO  
     200                     IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out) - ssh(ji,jj,Kbb_a) 
     201 
     202                     CALL remap_linear(tabin(1:N_in,1:jptra), z_in(1:N_in), tabres(ji,jj,1:N_out,1:jptra), & 
     203                                         &   z_out(1:N_out), N_in, N_out, jptra) 
     204                  END DO 
     205               END DO 
     206            ENDIF 
     207 
     208            DO jj=j1,j2 
     209               DO ji=i1,i2 
     210                  DO jk=1,jpkm1 
     211                     trbdiff(ji,jj,jk,1:jptra) = (tr(ji,jj,jk,1:jptra,Kbb_a) - tabres(ji,jj,jk,1:jptra))*tmask(ji,jj,jk) 
     212                  END DO 
     213               END DO 
     214            END DO 
     215 
     216         END IF 
     217 
     218         DO jn = 1, jptra             
     219            DO jk = 1, jpkm1 
     220               ztu(i1-1:i2,j1-1:j2,jk) = 0._wp 
     221               DO jj = j1,j2 
     222                  DO ji = i1,i2-1 
     223                     zabe1 = rn_sponge_tra * r1_Dt * umask(ji,jj,jk) * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) 
     224                     ztu(ji,jj,jk) = zabe1 * fspu(ji,jj) * ( trbdiff(ji+1,jj  ,jk,jn) - trbdiff(ji,jj,jk,jn) )  
     225                  END DO 
     226               END DO 
     227               ztv(i1-1:i2,j1-1:j2,jk) = 0._wp 
     228               DO ji = i1,i2 
     229                  DO jj = j1,j2-1 
     230                     zabe2 = rn_sponge_tra * r1_Dt * vmask(ji,jj,jk) * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a) 
     231                     ztv(ji,jj,jk) = zabe2 * fspv(ji,jj) * ( trbdiff(ji  ,jj+1,jk,jn) - trbdiff(ji,jj,jk,jn) ) 
     232                  END DO 
     233               END DO 
     234               ! 
     235               IF( ln_zps ) THEN      ! set gradient at partial step level 
     236                  DO jj = j1,j2 
     237                     DO ji = i1,i2 
     238                        ! last level 
     239                        iku = mbku(ji,jj) 
     240                        ikv = mbkv(ji,jj) 
     241                        IF( iku == jk )   ztu(ji,jj,jk) = 0._wp 
     242                        IF( ikv == jk )   ztv(ji,jj,jk) = 0._wp 
     243                     END DO 
     244                  END DO 
    121245               ENDIF 
    122             ENDDO 
    123          ENDDO 
    124 # endif 
    125  
    126          DO jj=j1,j2 
    127             DO ji=i1,i2 
    128                DO jk=1,jpkm1 
    129 # if defined key_vertical 
    130                   trbdiff(ji,jj,jk,1:jptra) = tr(ji,jj,jk,1:jptra,Kbb_a) - tabres_child(ji,jj,jk,1:jptra) 
    131 # else 
    132                   trbdiff(ji,jj,jk,1:jptra) = tr(ji,jj,jk,1:jptra,Kbb_a) - tabres(ji,jj,jk,1:jptra) 
    133 # endif 
    134                ENDDO 
    135             ENDDO 
    136          ENDDO 
    137  
    138          !* set relaxation time scale 
    139          IF( l_1st_euler .AND. lk_agrif_fstep ) THEN   ;   ztrelax =   rn_trelax_tra  / (        rn_Dt ) 
    140          ELSE                                          ;   ztrelax =   rn_trelax_tra  / (2._wp * rn_Dt ) 
    141          ENDIF 
    142  
    143          DO jn = 1, jptra 
     246            END DO 
     247            ! 
     248! JC: there is something wrong with the Laplacian in corners 
    144249            DO jk = 1, jpkm1 
    145                DO jj = j1,j2-1 
    146                   DO ji = i1,i2-1 
    147                      zabe1 = rn_sponge_tra * fspu(ji,jj) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * umask(ji,jj,jk) 
    148                      zabe2 = rn_sponge_tra * fspv(ji,jj) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) 
    149                      ztu(ji,jj) = zabe1 * ( trbdiff(ji+1,jj  ,jk,jn) - trbdiff(ji,jj,jk,jn) ) 
    150                      ztv(ji,jj) = zabe2 * ( trbdiff(ji  ,jj+1,jk,jn) - trbdiff(ji,jj,jk,jn) ) 
    151                   END DO 
    152                END DO 
    153                ! 
    154                DO jj = j1+1,j2-1 
    155                   DO ji = i1+1,i2-1 
    156                      IF( .NOT. tabspongedone_trn(ji,jj) ) THEN  
    157                         tr(ji,jj,jk,jn,Krhs_a) = tr(ji,jj,jk,jn,Krhs_a) + (  ztu(ji,jj) - ztu(ji-1,jj  )     & 
    158                            &                                   + ztv(ji,jj) - ztv(ji  ,jj-1)  )  & 
    159                            &                                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a)  & 
    160                            &                                - ztrelax * fspt(ji,jj) * trbdiff(ji,jj,jk,jn) 
     250               DO jj = j1,j2 
     251                  DO ji = i1,i2 
     252                     IF (.NOT. tabspongedone_trn(ji,jj)) THEN  
     253                        zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a) 
     254                        ! horizontal diffusive trends 
     255                        ztra = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   &  
     256                          &           + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) & 
     257                          &   - rn_trelax_tra * r1_Dt * fspt(ji,jj) * trbdiff(ji,jj,jk,jn) 
     258                        ! add it to the general tracer trends 
     259                        tr(ji,jj,jk,jn,Krhs_a) = tr(ji,jj,jk,jn,Krhs_a) + ztra 
    161260                     ENDIF 
    162261                  END DO 
    163262               END DO 
     263 
    164264            END DO 
    165265            ! 
    166266         END DO 
    167267         ! 
    168          tabspongedone_trn(i1+1:i2-1,j1+1:j2-1) = .TRUE. 
     268         tabspongedone_trn(i1:i2,j1:j2) = .TRUE. 
     269         ! 
    169270      ENDIF 
    170       !                  
     271      ! 
    171272   END SUBROUTINE interptrn_sponge 
    172273 

  • NEMO/trunk/src/NST/agrif_top_update.F90

    r12489 r14086  
    4040      IF (Agrif_Root()) RETURN  
    4141      ! 
    42       Agrif_UseSpecialValueInUpdate = .TRUE. 
     42      l_vremap                      = ln_vert_remap 
     43      Agrif_UseSpecialValueInUpdate = .NOT.l_vremap 
    4344      Agrif_SpecialValueFineGrid    = 0._wp 
     45 
    4446      !  
    4547# if ! defined DECAL_FEEDBACK 
    … …  
    5254      ! 
    5355      Agrif_UseSpecialValueInUpdate = .FALSE. 
     56      l_vremap                      = .FALSE. 
    5457      ! 
    5558   END SUBROUTINE Agrif_Update_Trc 
    5659 
    57 #ifdef key_vertical 
    5860   SUBROUTINE updateTRC( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    59       !!--------------------------------------------- 
    60       !!           *** ROUTINE updateT *** 
    61       !!--------------------------------------------- 
     61 
    6262      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2 
    6363      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    … …  
    7272      REAL(wp) :: tabin(k1:k2,1:jptra) 
    7373      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,1:jptra) :: tabres_child 
    74       !!--------------------------------------------- 
    75       ! 
     74 
    7675      IF (before) THEN 
    77          AGRIF_SpecialValue = -999._wp 
    78          DO jn = n1,n2-1 
     76         IF ( l_vremap ) THEN 
     77            DO jn = n1,n2-1 
     78               DO jk=k1,k2 
     79                  DO jj=j1,j2 
     80                     DO ji=i1,i2 
     81                        tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) 
     82                     END DO 
     83                  END DO 
     84               END DO 
     85            END DO 
    7986            DO jk=k1,k2 
    8087               DO jj=j1,j2 
    8188                  DO ji=i1,i2 
    82                      tabres(ji,jj,jk,jn) = (tr(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) ) & 
    83                                            * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp 
    84                   END DO 
    85                END DO 
    86             END DO 
    87          END DO 
    88          DO jk=k1,k2 
     89                     tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 
     90                  END DO 
     91               END DO 
     92            END DO 
     93         ELSE 
     94            DO jn = 1,jptra 
     95               DO jk=k1,k2 
     96                  DO jj=j1,j2 
     97                     DO ji=i1,i2 
     98                        tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kmm_a)  * e3t(ji,jj,jk,Kmm_a) / e3t_0(ji,jj,jk) 
     99                     END DO 
     100                  END DO 
     101               END DO 
     102            END DO 
     103 
     104         ENDIF 
     105      ELSE 
     106         IF ( l_vremap ) THEN 
     107            tabres_child(:,:,:,:) = 0._wp 
     108            AGRIF_SpecialValue = 0._wp 
    89109            DO jj=j1,j2 
    90110               DO ji=i1,i2 
    91                   tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) & 
    92                                            + (tmask(ji,jj,jk)-1)*999._wp 
    93                END DO 
    94             END DO 
    95          END DO 
    96       ELSE 
    97          tabres_child(:,:,:,:) = 0. 
    98          AGRIF_SpecialValue = 0._wp 
    99          DO jj=j1,j2 
    100             DO ji=i1,i2 
    101                N_in = 0 
    102                DO jk=k1,k2 !k2 = jpk of child grid 
    103                   IF (tabres(ji,jj,jk,n2) == 0  ) EXIT 
    104                   N_in = N_in + 1 
    105                   tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 
    106                   h_in(N_in) = tabres(ji,jj,jk,n2) 
     111                  N_in = 0 
     112                  DO jk=k1,k2 !k2 = jpk of child grid 
     113                     IF (tabres(ji,jj,jk,n2) <= 1.e-6_wp  ) EXIT 
     114                     N_in = N_in + 1 
     115                     tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 
     116                     h_in(N_in) = tabres(ji,jj,jk,n2) 
     117                  ENDDO 
     118                  N_out = 0 
     119                  DO jk=1,jpk ! jpk of parent grid 
     120                     IF (tmask(ji,jj,jk) == 0 ) EXIT ! TODO: Will not work with ISF 
     121                     N_out = N_out + 1 
     122                     h_out(N_out) = e3t(ji,jj,jk,Kmm_a) 
     123                  ENDDO 
     124                  IF (N_in*N_ou