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 13937 for NEMO/branches – NEMO

Changeset 13937 for NEMO/branches


Ignore:
Timestamp:
2020-12-01T13:44:17+01:00 (3 years ago)
Author:
jchanut
Message:

#2222, added:

  • Divergence conserving interpolation for transport (Balsara's scheme)
  • Possibility to shift barotropic interface inside the domain (currently disabled - needed to match volume within sponge)
  • PPR library as default vertical interpolation scheme
  • Trilinear interpolation over partial bottom cells
Location:
NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_oce.F90

    r13565 r13937  
    3131   ! 
    3232   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 
    3334 
    3435   LOGICAL , PUBLIC :: spongedoneT = .FALSE.       !: tracer   sponge layer indicator 
     
    3637   LOGICAL , PUBLIC :: lk_agrif_fstep = .TRUE.     !: if true: first step 
    3738   LOGICAL , PUBLIC :: lk_agrif_debug = .FALSE.    !: if true: print debugging info 
    38  
     39   LOGICAL , PUBLIC :: lk_tint2d_notinterp = .FALSE. !: if true, no time interp 
    3940   LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_tsn 
    4041# if defined key_top 
     
    4748   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fspu, fspv !: sponge arrays 
    4849   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 !:   "       "     "   " 
    4952 
    5053   ! Barotropic arrays used to store open boundary data during time-splitting loop: 
     
    5255   INTEGER , PUBLIC,              SAVE                 ::  Kbb_a, Kmm_a, Krhs_a   !: AGRIF module-specific copies of time-level indices 
    5356 
    54    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht0_parent, hu0_parent, hv0_parent 
    55    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 
     60 
    5661 
    5762   INTEGER, PUBLIC :: ts_interp_id, ts_update_id                              ! AGRIF profile for tracers interpolation and update 
    5863   INTEGER, PUBLIC :: un_interp_id, vn_interp_id                              ! AGRIF profiles for interpolations 
    5964   INTEGER, PUBLIC :: un_update_id, vn_update_id                              ! AGRIF profiles for udpates 
    60    INTEGER, PUBLIC :: ts_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) 
    6167   INTEGER, PUBLIC :: tsini_id, uini_id, vini_id, sshini_id                   ! AGRIF profile for initialization 
    6268# if defined key_top 
    6369   INTEGER, PUBLIC :: trn_id, trn_sponge_id 
    6470# endif   
    65    INTEGER, PUBLIC :: unb_id, vnb_id, ub2b_interp_id, vb2b_interp_id 
    66    INTEGER, PUBLIC :: ub2b_update_id, vb2b_update_id 
    67    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 
    6875   INTEGER, PUBLIC :: scales_t_id 
    6976   INTEGER, PUBLIC :: avt_id, avm_id, en_id                ! TKE related identificators 
    70    INTEGER, PUBLIC :: mbkt_id, ht0_id 
     77   INTEGER, PUBLIC :: mbkt_id, ht0_id, e3t0_interp_id 
    7178   INTEGER, PUBLIC :: glamt_id, gphit_id 
    7279   INTEGER, PUBLIC :: batupd_id 
     
    98105      ALLOCATE( fspu(jpi,jpj), fspv(jpi,jpj),          & 
    99106         &      fspt(jpi,jpj), fspf(jpi,jpj),               & 
     107         &      fspu_2d(jpi,jpj), fspv_2d(jpi,jpj),         & 
     108         &      fspt_2d(jpi,jpj), fspf_2d(jpi,jpj),         & 
    100109         &      tabspongedone_tsn(jpi,jpj),                 & 
    101110         &      utint_stage(jpi,jpj), vtint_stage(jpi,jpj), & 
  • NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_oce_interp.F90

    r13498 r13937  
    4545   PUBLIC   interpunb, interpvnb , interpub2b, interpvb2b 
    4646   PUBLIC   interpe3t, interpglamt, interpgphit 
    47    PUBLIC   interpht0, interpmbkt 
     47   PUBLIC   interpht0, interpmbkt, interpe3t0_vremap 
    4848   PUBLIC   agrif_istate_oce, agrif_istate_ssh   ! called by icestate.F90 and domvvl.F90 
    4949   PUBLIC   agrif_check_bat 
     
    200200      CALL Agrif_Bc_variable( un_interp_id, procname=interpun ) 
    201201      CALL Agrif_Bc_variable( vn_interp_id, procname=interpvn ) 
     202 
     203      IF( .NOT.ln_dynspg_ts ) THEN ! Get transports 
     204         ubdy(:,:) = 0._wp   ;   vbdy(:,:) = 0._wp  
     205         CALL Agrif_Bc_variable( unb_interp_id, procname=interpunb ) 
     206         CALL Agrif_Bc_variable( vnb_interp_id, procname=interpvnb ) 
     207      ENDIF 
     208 
    202209      use_sign_north = .FALSE. 
    203210      ! 
    204211      Agrif_UseSpecialValue = .FALSE. 
    205212      l_vremap              = .FALSE. 
     213      ! 
     214      ! Ensure below that vertically integrated transports match 
     215      ! either transports out of time splitting procedure (ln_dynspg_ts=.TRUE.) 
     216      ! or parent grid transports (ln_dynspg_ts=.FALSE.) 
    206217      ! 
    207218      ! --- West --- ! 
    208219      IF( lk_west ) THEN 
    209220         ibdy1 = nn_hls + 2                  ! halo + land + 1 
    210          ibdy2 = nn_hls + 1 + nbghostcells   ! halo + land + nbghostcells 
     221         ibdy2 = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhox()   ! halo + land + nbghostcells 
    211222         ! 
    212223         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport 
    213224            DO ji = mi0(ibdy1), mi1(ibdy2) 
    214                uu_b(ji,:,Krhs_a) = 0._wp 
    215                DO jk = 1, jpkm1 
    216                   DO jj = 1, jpj 
    217                      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) 
    218                   END DO 
    219                END DO 
    220225               DO jj = 1, jpj 
    221                   uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) * r1_hu(ji,jj,Krhs_a) 
     226                  uu_b(ji,jj,Krhs_a) = ubdy(ji,jj) * r1_hu(ji,jj,Krhs_a) 
     227                  vv_b(ji,jj,Krhs_a) = vbdy(ji,jj) * r1_hv(ji,jj,Krhs_a) 
    222228               END DO 
    223229            END DO 
     
    225231         ! 
    226232         DO ji = mi0(ibdy1), mi1(ibdy2) 
    227             zub(ji,:) = 0._wp    ! Correct transport 
     233            zub(ji,:) = 0._wp   
    228234            DO jk = 1, jpkm1 
    229235               DO jj = 1, jpj 
     
    241247         END DO 
    242248         !    
    243          IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate 
    244             DO ji = mi0(ibdy1), mi1(ibdy2) 
    245                zvb(ji,:) = 0._wp 
    246                DO jk = 1, jpkm1 
    247                   DO jj = 1, jpj 
    248                      zvb(ji,jj) = zvb(ji,jj) + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 
    249                   END DO 
    250                END DO 
     249         DO ji = mi0(ibdy1), mi1(ibdy2) 
     250            zvb(ji,:) = 0._wp 
     251            DO jk = 1, jpkm1 
    251252               DO jj = 1, jpj 
    252                   zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a) 
    253                END DO 
    254                DO jk = 1, jpkm1 
    255                   DO jj = 1, jpj 
    256                      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) 
    257                   END DO 
    258                END DO 
    259             END DO 
    260          ENDIF 
     253                  zvb(ji,jj) = zvb(ji,jj) + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 
     254               END DO 
     255            END DO 
     256            DO jj = 1, jpj 
     257               zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a) 
     258            END DO 
     259            DO jk = 1, jpkm1 
     260               DO jj = 1, jpj 
     261                  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) 
     262               END DO 
     263            END DO 
     264         END DO 
    261265         ! 
    262266      ENDIF 
     
    264268      ! --- East --- ! 
    265269      IF( lk_east) THEN 
    266          ibdy1 = jpiglo - ( nn_hls + nbghostcells + 1)   ! halo + land + nbghostcells 
    267          ibdy2 = jpiglo - ( nn_hls + 2 )                 ! halo + land + 1 
    268          ! 
    269          IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport 
     270         ibdy1 = jpiglo - ( nn_hls + nbghostcells + 1) - nn_shift_bar*Agrif_Rhox()     
     271         ibdy2 = jpiglo - ( nn_hls + 2 )                  
     272         ! 
     273         IF( .NOT.ln_dynspg_ts ) THEN  
    270274            DO ji = mi0(ibdy1), mi1(ibdy2) 
    271275               uu_b(ji,:,Krhs_a) = 0._wp 
     
    276280               END DO 
    277281               DO jj = 1, jpj 
    278                   uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) * r1_hu(ji,jj,Krhs_a) 
     282                  uu_b(ji,jj,Krhs_a) = ubdy(ji,jj) * r1_hu(ji,jj,Krhs_a) 
    279283               END DO 
    280284            END DO 
     
    282286         ! 
    283287         DO ji = mi0(ibdy1), mi1(ibdy2) 
    284             zub(ji,:) = 0._wp    ! Correct transport 
     288            zub(ji,:) = 0._wp    
    285289            DO jk = 1, jpkm1 
    286290               DO jj = 1, jpj 
     
    298302         END DO 
    299303         ! 
    300          IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate 
    301             ibdy1 = jpiglo - ( nn_hls + nbghostcells )   ! halo + land + nbghostcells - 1 
    302             ibdy2 = jpiglo - ( nn_hls + 1 )              ! halo + land + 1            - 1 
     304         ibdy1 = jpiglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhox()  
     305         ibdy2 = jpiglo - ( nn_hls + 1 )              ! 
     306         IF( .NOT.ln_dynspg_ts ) THEN  
    303307            DO ji = mi0(ibdy1), mi1(ibdy2) 
    304                zvb(ji,:) = 0._wp 
     308               vv_b(ji,:,Krhs_a) = 0._wp 
    305309               DO jk = 1, jpkm1 
    306310                  DO jj = 1, jpj 
     311                     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) 
     312                  END DO 
     313               END DO 
     314               DO jj = 1, jpj 
     315                  vv_b(ji,jj,Krhs_a) = vbdy(ji,jj) * r1_hv(ji,jj,Krhs_a) 
     316               END DO 
     317            END DO 
     318         ENDIF 
     319 
     320         DO ji = mi0(ibdy1), mi1(ibdy2) 
     321            zvb(ji,:) = 0._wp 
     322            DO jk = 1, jpkm1 
     323               DO jj = 1, jpj 
    307324                     zvb(ji,jj) = zvb(ji,jj) + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 
    308                   END DO 
    309                END DO 
     325               END DO 
     326            END DO 
     327            DO jj = 1, jpj 
     328               zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a) 
     329            END DO 
     330            DO jk = 1, jpkm1 
    310331               DO jj = 1, jpj 
    311                   zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a) 
    312                END DO 
    313                DO jk = 1, jpkm1 
    314                   DO jj = 1, jpj 
    315332                     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) 
    316                   END DO 
    317                END DO 
    318             END DO 
    319          ENDIF 
     333               END DO 
     334            END DO 
     335         END DO 
    320336         ! 
    321337      ENDIF 
     
    323339      ! --- South --- ! 
    324340      IF( lk_south ) THEN 
    325          jbdy1 = nn_hls + 2                  ! halo + land + 1 
    326          jbdy2 = nn_hls + 1 + nbghostcells   ! halo + land + nbghostcells 
    327          ! 
    328          IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport 
     341         jbdy1 = nn_hls + 2                  
     342         jbdy2 = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhoy()    
     343         ! 
     344         IF( .NOT.ln_dynspg_ts ) THEN 
    329345            DO jj = mj0(jbdy1), mj1(jbdy2) 
    330346               vv_b(:,jj,Krhs_a) = 0._wp 
     
    335351               END DO 
    336352               DO ji=1,jpi 
    337                   vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) * r1_hv(ji,jj,Krhs_a)      
     353                  vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) * r1_hv(ji,jj,Krhs_a)  
     354                  uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) * r1_hu(ji,jj,Krhs_a) 
    338355               END DO 
    339356            END DO 
     
    341358         ! 
    342359         DO jj = mj0(jbdy1), mj1(jbdy2) 
    343             zvb(:,jj) = 0._wp    ! Correct transport 
     360            zvb(:,jj) = 0._wp 
    344361            DO jk=1,jpkm1 
    345362               DO ji=1,jpi 
     
    357374         END DO 
    358375         ! 
    359          IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate 
    360             DO jj = mj0(jbdy1), mj1(jbdy2) 
    361                zub(:,jj) = 0._wp 
    362                DO jk = 1, jpkm1 
    363                   DO ji = 1, jpi 
    364                      zub(ji,jj) = zub(ji,jj) + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 
    365                   END DO 
    366                END DO 
     376         DO jj = mj0(jbdy1), mj1(jbdy2) 
     377            zub(:,jj) = 0._wp 
     378            DO jk = 1, jpkm1 
    367379               DO ji = 1, jpi 
    368                   zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 
    369                END DO 
    370                DO jk = 1, jpkm1 
    371                   DO ji = 1, jpi 
    372                      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) 
    373                   END DO 
    374                END DO 
    375             END DO 
    376          ENDIF 
     380                  zub(ji,jj) = zub(ji,jj) + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 
     381               END DO 
     382            END DO 
     383            DO ji = 1, jpi 
     384               zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 
     385            END DO 
     386            DO jk = 1, jpkm1 
     387               DO ji = 1, jpi 
     388                  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) 
     389               END DO 
     390            END DO 
     391         END DO 
    377392         ! 
    378393      ENDIF 
     
    380395      ! --- North --- ! 
    381396      IF( lk_north ) THEN 
    382          jbdy1 = jpjglo - ( nn_hls + nbghostcells + 1)   ! halo + land + nbghostcells 
    383          jbdy2 = jpjglo - ( nn_hls + 2 )                 ! halo + land + 1 
    384          ! 
    385          IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport 
     397         jbdy1 = jpjglo - ( nn_hls + nbghostcells + 1) - nn_shift_bar*Agrif_Rhoy()  
     398         jbdy2 = jpjglo - ( nn_hls + 2 ) 
     399         ! 
     400         IF( .NOT.ln_dynspg_ts ) THEN 
    386401            DO jj = mj0(jbdy1), mj1(jbdy2) 
    387402               vv_b(:,jj,Krhs_a) = 0._wp 
     
    392407               END DO 
    393408               DO ji=1,jpi 
    394                   vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) * r1_hv(ji,jj,Krhs_a) 
     409                  vv_b(ji,jj,Krhs_a) = vbdy(ji,jj) * r1_hv(ji,jj,Krhs_a) 
    395410               END DO 
    396411            END DO 
     
    398413         ! 
    399414         DO jj = mj0(jbdy1), mj1(jbdy2) 
    400             zvb(:,jj) = 0._wp    ! Correct transport 
     415            zvb(:,jj) = 0._wp  
    401416            DO jk=1,jpkm1 
    402417               DO ji=1,jpi 
     
    414429         END DO 
    415430         ! 
    416          IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate 
    417             jbdy1 = jpjglo - ( nn_hls + nbghostcells )   ! halo + land + nbghostcells - 1 
    418             jbdy2 = jpjglo - ( nn_hls + 1 )              ! halo + land + 1            - 1 
     431         jbdy1 = jpjglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhoy()   
     432         jbdy2 = jpjglo - ( nn_hls + 1 ) 
     433         IF( .NOT.ln_dynspg_ts ) THEN 
    419434            DO jj = mj0(jbdy1), mj1(jbdy2) 
    420                zub(:,jj) = 0._wp 
     435               uu_b(:,jj,Krhs_a) = 0._wp 
    421436               DO jk = 1, jpkm1 
    422437                  DO ji = 1, jpi 
    423                      zub(ji,jj) = zub(ji,jj) + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 
    424                   END DO 
    425                END DO 
     438                     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) 
     439                  END DO 
     440               END DO 
     441               DO ji=1,jpi 
     442                  uu_b(ji,jj,Krhs_a) = ubdy(ji,jj) * r1_hu(ji,jj,Krhs_a) 
     443               END DO 
     444            END DO 
     445         ENDIF 
     446 
     447         DO jj = mj0(jbdy1), mj1(jbdy2) 
     448            zub(:,jj) = 0._wp 
     449            DO jk = 1, jpkm1 
    426450               DO ji = 1, jpi 
    427                   zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 
    428                END DO 
    429                DO jk = 1, jpkm1 
    430                   DO ji = 1, jpi 
     451                  zub(ji,jj) = zub(ji,jj) + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 
     452               END DO 
     453            END DO 
     454            DO ji = 1, jpi 
     455               zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 
     456            END DO 
     457            DO jk = 1, jpkm1 
     458               DO ji = 1, jpi 
    431459                     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) 
    432                   END DO 
    433                END DO 
    434             END DO 
    435          ENDIF 
     460               END DO 
     461            END DO 
     462         END DO 
    436463         ! 
    437464      ENDIF 
     
    455482      IF( lk_west ) THEN 
    456483         istart = nn_hls + 2                              ! halo + land + 1 
    457          iend   = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
     484         iend   = nn_hls + 1 + nbghostcells  + nn_shift_bar*Agrif_Rhox()              ! halo + land + nbghostcells 
    458485         DO ji = mi0(istart), mi1(iend) 
    459486            DO jj=1,jpj 
     
    466493      !--- East ---! 
    467494      IF( lk_east ) THEN 
    468          istart = jpiglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    469          iend   = jpiglo - ( nn_hls + 1 )                 ! halo + land + 1            - 1 
     495         istart = jpiglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhox()  
     496         iend   = jpiglo - ( nn_hls + 1 )                 
    470497         DO ji = mi0(istart), mi1(iend) 
    471498 
     
    474501            END DO 
    475502         END DO 
    476          istart = jpiglo - ( nn_hls + nbghostcells + 1)   ! halo + land + nbghostcells 
    477          iend   = jpiglo - ( nn_hls + 2 )                 ! halo + land + 1 
     503         istart = jpiglo - ( nn_hls + nbghostcells + 1) - nn_shift_bar*Agrif_Rhox()  
     504         iend   = jpiglo - ( nn_hls + 2 )                 
    478505         DO ji = mi0(istart), mi1(iend) 
    479506            DO jj=1,jpj 
     
    485512      !--- South ---! 
    486513      IF( lk_south ) THEN 
    487          jstart = nn_hls + 2                              ! halo + land + 1 
    488          jend   = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
     514         jstart = nn_hls + 2                               
     515         jend   = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhoy()            
    489516         DO jj = mj0(jstart), mj1(jend) 
    490517 
     
    498525      !--- North ---! 
    499526      IF( lk_north ) THEN 
    500          jstart = jpjglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    501          jend   = jpjglo - ( nn_hls + 1 )                 ! halo + land + 1            - 1 
     527         jstart = jpjglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhoy()      
     528         jend   = jpjglo - ( nn_hls + 1 )                 
    502529         DO jj = mj0(jstart), mj1(jend) 
    503530            DO ji=1,jpi 
     
    505532            END DO 
    506533         END DO 
    507          jstart = jpjglo - ( nn_hls + nbghostcells + 1)   ! halo + land + nbghostcells 
    508          jend   = jpjglo - ( nn_hls + 2 )                 ! halo + land + 1 
     534         jstart = jpjglo - ( nn_hls + nbghostcells + 1) - nn_shift_bar*Agrif_Rhoy()  
     535         jend   = jpjglo - ( nn_hls + 2 )                 
    509536         DO jj = mj0(jstart), mj1(jend) 
    510537            DO ji=1,jpi 
     
    532559      !--- West ---! 
    533560      IF( lk_west ) THEN 
    534          istart = nn_hls + 2                              ! halo + land + 1 
    535          iend   = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
     561         istart = nn_hls + 2                               
     562         iend   = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhox()  
    536563         DO ji = mi0(istart), mi1(iend) 
    537564            DO jj=1,jpj 
     
    544571      !--- East ---! 
    545572      IF( lk_east ) THEN 
    546          istart = jpiglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    547          iend   = jpiglo - ( nn_hls + 1 )                 ! halo + land + 1            - 1 
     573         istart = jpiglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhox() 
     574         iend   = jpiglo - ( nn_hls + 1 )                  
    548575         DO ji = mi0(istart), mi1(iend) 
    549576            DO jj=1,jpj 
     
    551578            END DO 
    552579         END DO 
    553          istart = jpiglo - ( nn_hls + nbghostcells + 1)   ! halo + land + nbghostcells 
    554          iend   = jpiglo - ( nn_hls + 2 )                 ! halo + land + 1 
     580         istart = jpiglo - ( nn_hls + nbghostcells + 1) - nn_shift_bar*Agrif_Rhox()  
     581         iend   = jpiglo - ( nn_hls + 2 )                  
    555582         DO ji = mi0(istart), mi1(iend) 
    556583            DO jj=1,jpj 
     
    562589      !--- South ---! 
    563590      IF( lk_south ) THEN 
    564          jstart = nn_hls + 2                              ! halo + land + 1 
    565          jend   = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
     591         jstart = nn_hls + 2                               
     592         jend   = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhoy()  
    566593         DO jj = mj0(jstart), mj1(jend) 
    567594            DO ji=1,jpi 
     
    574601      !--- North ---! 
    575602      IF( lk_north ) THEN 
    576          jstart = jpjglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    577          jend   = jpjglo - ( nn_hls + 1 )                 ! halo + land + 1            - 1 
     603         jstart = jpjglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhoy()  
     604         jend   = jpjglo - ( nn_hls + 1 )                 
    578605         DO jj = mj0(jstart), mj1(jend) 
    579606            DO ji=1,jpi 
     
    581608            END DO 
    582609         END DO 
    583          jstart = jpjglo - ( nn_hls + nbghostcells + 1)   ! halo + land + nbghostcells 
    584          jend   = jpjglo - ( nn_hls + 2 )                 ! halo + land + 1 
     610         jstart = jpjglo - ( nn_hls + nbghostcells + 1) - nn_shift_bar*Agrif_Rhoy()  
     611         jend   = jpjglo - ( nn_hls + 2 )                
    585612         DO jj = mj0(jstart), mj1(jend) 
    586613            DO ji=1,jpi 
     
    599626      INTEGER, INTENT(in) ::   kt 
    600627      !! 
    601       INTEGER :: ji, jj 
    602628      LOGICAL :: ll_int_cons 
    603629      !!----------------------------------------------------------------------   
     
    623649      ! 
    624650      IF( ll_int_cons ) THEN  ! Conservative interpolation 
    625          ! order matters here !!!!!! 
    626          CALL Agrif_Bc_variable( ub2b_interp_id, calledweight=1._wp, procname=interpub2b ) ! Time integrated 
    627          CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b )  
    628          ! 
    629          bdy_tinterp = 1 
    630          CALL Agrif_Bc_variable( unb_id        , calledweight=1._wp, procname=interpunb  ) ! After 
    631          CALL Agrif_Bc_variable( vnb_id        , calledweight=1._wp, procname=interpvnb  )   
    632          ! 
    633          bdy_tinterp = 2 
    634          CALL Agrif_Bc_variable( unb_id        , calledweight=0._wp, procname=interpunb  ) ! Before 
    635          CALL Agrif_Bc_variable( vnb_id        , calledweight=0._wp, procname=interpvnb  )    
     651         IF ( lk_tint2d_notinterp ) THEN 
     652            CALL Agrif_Bc_variable( ub2b_interp_id, calledweight=1._wp, procname=interpub2b_const ) 
     653            CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b_const )  
     654            ! Divergence conserving correction terms: 
     655            CALL Agrif_Bc_variable(    ub2b_cor_id, calledweight=1._wp, procname=        ub2b_cor ) 
     656            CALL Agrif_Bc_variable(    vb2b_cor_id, calledweight=1._wp, procname=        vb2b_cor ) 
     657         ELSE 
     658            ! order matters here !!!!!! 
     659            CALL Agrif_Bc_variable( ub2b_interp_id, calledweight=1._wp, procname=interpub2b ) ! Time integrated 
     660            CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b )  
     661            ! 
     662            bdy_tinterp = 1 
     663            CALL Agrif_Bc_variable( unb_interp_id , calledweight=1._wp, procname=interpunb  ) ! After 
     664            CALL Agrif_Bc_variable( vnb_interp_id , calledweight=1._wp, procname=interpvnb  )   
     665            ! 
     666            bdy_tinterp = 2 
     667            CALL Agrif_Bc_variable( unb_interp_id , calledweight=0._wp, procname=interpunb  ) ! Before 
     668            CALL Agrif_Bc_variable( vnb_interp_id , calledweight=0._wp, procname=interpvnb  )    
     669         ENDIF 
    636670      ELSE ! Linear interpolation 
    637671         ! 
    638672         ubdy(:,:) = 0._wp   ;   vbdy(:,:) = 0._wp  
    639          CALL Agrif_Bc_variable( unb_id, procname=interpunb ) 
    640          CALL Agrif_Bc_variable( vnb_id, procname=interpvnb ) 
     673         CALL Agrif_Bc_variable( unb_interp_id, procname=interpunb ) 
     674         CALL Agrif_Bc_variable( vnb_interp_id, procname=interpvnb ) 
    641675      ENDIF 
    642676      Agrif_UseSpecialValue = .FALSE. 
     
    667701      ! --- West --- ! 
    668702      IF(lk_west) THEN 
    669          istart = nn_hls + 2                              ! halo + land + 1 
    670          iend   = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
     703         istart = nn_hls + 2                                                          ! halo + land + 1 
     704         iend   = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhox()               ! halo + land + nbghostcells 
    671705         DO ji = mi0(istart), mi1(iend) 
    672706            DO jj = 1, jpj 
     
    678712      ! --- East --- ! 
    679713      IF(lk_east) THEN 
    680          istart = jpiglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    681          iend   = jpiglo - ( nn_hls + 1 )                 ! halo + land + 1            - 1 
     714         istart = jpiglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhox()       ! halo + land + nbghostcells - 1 
     715         iend   = jpiglo - ( nn_hls + 1 )                                              ! halo + land + 1            - 1 
    682716         DO ji = mi0(istart), mi1(iend) 
    683717            DO jj = 1, jpj 
     
    689723      ! --- South --- ! 
    690724      IF(lk_south) THEN 
    691          jstart = nn_hls + 2                              ! halo + land + 1 
    692          jend   = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
     725         jstart = nn_hls + 2                                                          ! halo + land + 1 
     726         jend   = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhoy()               ! halo + land + nbghostcells 
    693727         DO jj = mj0(jstart), mj1(jend) 
    694728            DO ji = 1, jpi 
     
    700734      ! --- North --- ! 
    701735      IF(lk_north) THEN 
    702          jstart = jpjglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    703          jend   = jpjglo - ( nn_hls + 1 )                 ! halo + land + 1            - 1 
     736         jstart = jpjglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhoy()     ! halo + land + nbghostcells - 1 
     737         jend   = jpjglo - ( nn_hls + 1 )                                            ! halo + land + 1            - 1 
    704738         DO jj = mj0(jstart), mj1(jend) 
    705739            DO ji = 1, jpi 
     
    726760      ! --- West --- ! 
    727761      IF(lk_west) THEN 
    728          istart = nn_hls + 2                              ! halo + land + 1 
    729          iend   = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
     762         istart = nn_hls + 2                                                        ! halo + land + 1 
     763         iend   = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhox()             ! halo + land + nbghostcells 
    730764         DO ji = mi0(istart), mi1(iend) 
    731765            DO jj = 1, jpj 
     
    737771      ! --- East --- ! 
    738772      IF(lk_east) THEN 
    739          istart = jpiglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    740          iend   = jpiglo - ( nn_hls + 1 )                 ! halo + land + 1            - 1 
     773         istart = jpiglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhox()    ! halo + land + nbghostcells - 1 
     774         iend   = jpiglo - ( nn_hls + 1 )                                           ! halo + land + 1            - 1 
    741775         DO ji = mi0(istart), mi1(iend) 
    742776            DO jj = 1, jpj 
     
    748782      ! --- South --- ! 
    749783      IF(lk_south) THEN 
    750          jstart = nn_hls + 2                              ! halo + land + 1 
    751          jend   = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
     784         jstart = nn_hls + 2                                                        ! halo + land + 1 
     785         jend   = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhoy()             ! halo + land + nbghostcells 
    752786         DO jj = mj0(jstart), mj1(jend) 
    753787            DO ji = 1, jpi 
     
    759793      ! --- North --- ! 
    760794      IF(lk_north) THEN 
    761          jstart = jpjglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    762          jend   = jpjglo - ( nn_hls + 1 )                 ! halo + land + 1            - 1 
     795         jstart = jpjglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhoy()    ! halo + land + nbghostcells - 1 
     796         jend   = jpjglo - ( nn_hls + 1 )                                           ! halo + land + 1            - 1 
    763797         DO jj = mj0(jstart), mj1(jend) 
    764798            DO ji = 1, jpi 
     
    807841      INTEGER  :: item 
    808842      ! vertical interpolation: 
    809       REAL(wp) :: zhtot 
    810       REAL(wp), DIMENSION(k1:k2,1:jpts) :: tabin 
    811       REAL(wp), DIMENSION(k1:k2) :: h_in, z_in 
     843      REAL(wp) :: zhtot, zwgt 
     844      REAL(wp), DIMENSION(k1:k2,1:jpts) :: tabin, tabin_i 
     845      REAL(wp), DIMENSION(k1:k2) :: z_in, h_in_i, z_in_i 
    812846      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out 
    813847      !!---------------------------------------------------------------------- 
     
    828862         END DO 
    829863 
    830          IF( l_vremap .OR. l_ini_child ) THEN 
    831             ! Interpolate thicknesses 
     864         IF( l_vremap .OR. l_ini_child .OR. ln_zps ) THEN 
     865 
     866            ! Fill cell depths (i.e. gdept) to be interpolated 
    832867            ! Warning: these are masked, hence extrapolated prior interpolation. 
    833             DO jk=k1,k2 
    834                DO jj=j1,j2 
    835                   DO ji=i1,i2 
    836                       ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 
    837  
    838                   END DO 
    839                END DO 
    840             END DO 
    841  
    842             ! Extrapolate thicknesses in partial bottom cells: 
    843             ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
    844             IF (ln_zps) THEN 
    845                DO jj=j1,j2 
    846                   DO ji=i1,i2 
    847                       jk = mbkt(ji,jj) 
    848                       ptab(ji,jj,jk,jpts+1) = 0._wp 
    849                   END DO 
    850                END DO            
    851             END IF 
    852          
     868            DO jj=j1,j2 
     869               DO ji=i1,i2 
     870                  ptab(ji,jj,k1,jpts+1) = 0.5_wp * tmask(ji,jj,k1) * e3t(ji,jj,k1,Kmm_a) 
     871                  DO jk=k1+1,k2 
     872                     ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * & 
     873                        & ( ptab(ji,jj,jk-1,jpts+1) + 0.5_wp * (e3t(ji,jj,jk-1,Kmm_a)+e3t(ji,jj,jk,Kmm_a)) ) 
     874                  END DO 
     875               END DO 
     876            END DO 
     877          
    853878            ! Save ssh at last level: 
    854879            IF (.NOT.ln_linssh) THEN 
    855880               ptab(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1)  
    856             ELSE 
    857                ptab(i1:i2,j1:j2,k2,jpts+1) = 0._wp 
    858881            END IF       
    859882         ENDIF 
     
    866889         IF( l_vremap .OR. l_ini_child ) THEN 
    867890            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,n2) = 0._wp  
    868                 
    869891            DO jj=j1,j2 
    870892               DO ji=i1,i2 
    871                   ts(ji,jj,:,:,Krhs_a) = 0.                   
     893                  ts(ji,jj,:,:,Krhs_a) = 0.   
     894                  ! 
     895                  ! Build vertical grids: 
    872896                  N_in = mbkt_parent(ji,jj) 
    873                   zhtot = 0._wp 
    874                   DO jk=1,N_in !k2 = jpk of parent grid 
    875                      IF (jk==N_in) THEN 
    876                         h_in(jk) = ht0_parent(ji,jj) + ptab(ji,jj,k2,n2) - zhtot 
    877                      ELSE 
    878                         h_in(jk) = ptab(ji,jj,jk,n2) 
    879                      ENDIF 
    880                      zhtot = zhtot + h_in(jk) 
    881                      tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1) 
    882                   END DO 
    883                   z_in(1) = 0.5_wp * h_in(1) - zhtot + ht0_parent(ji,jj) 
     897                  ! Input grid (account for partial cells if any): 
     898                  DO jk=1,N_in 
     899                     z_in(jk) = ptab(ji,jj,jk,n2) - ptab(ji,jj,k2,n2) 
     900                     tabin(jk,1:jpts) = ptab(ji,jj,jk,1:jpts) 
     901                  END DO 
     902                   
     903                  ! Intermediate grid: 
     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) 
    884909                  DO jk=2,N_in 
    885                      z_in(jk) = z_in(jk-1) + 0.5_wp * (h_in(jk-1)+h_in(jk)) 
    886                   END DO 
    887  
    888                   N_out = 0 
    889                   DO jk=1,jpk ! jpk of child grid 
    890                      IF (tmask(ji,jj,jk) == 0._wp) EXIT  
    891                      N_out = N_out + 1 
     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 
     914                  ! Output (Child) grid: 
     915                  N_out = mbkt(ji,jj) 
     916                  DO jk=1,N_out 
    892917                     h_out(jk) = e3t(ji,jj,jk,Krhs_a) 
    893918                  END DO 
    894  
    895                   z_out(1) = 0.5_wp * h_out(1) - SUM(h_out(1:N_out)) + ht_0(ji,jj) 
     919                  z_out(1) = 0.5_wp * h_out(1) 
    896920                  DO jk=2,N_out 
    897                      z_out(jk) = z_out(jk-1) + 0.5_wp * (h_out(jk)+h_out(jk-1)) 
    898                   END DO 
     921                     z_out(jk) = z_out(jk-1) + 0.5_wp * ( h_out(jk)+h_out(jk-1) ) 
     922                  END DO 
     923                  IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out)  - ssh(ji,jj,Krhs_a) 
    899924 
    900925                  IF (N_in*N_out > 0) THEN 
    901926                     IF( l_ini_child ) THEN 
    902                         CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),          & 
     927                        CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),              & 
    903928                                      &   z_out(1:N_out),N_in,N_out,jpts)   
    904929                     ELSE  
    905                         CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),   & 
     930                        CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),tabin_i(1:N_in,1:jpts),                       & 
     931                                     &   z_in_i(1:N_in),N_in,N_in,jpts) 
     932                        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),   & 
    906933                                      &   h_out(1:N_out),N_in,N_out,jpts)   
    907934                     ENDIF 
     
    913940         ELSE 
    914941          
    915             DO jn=1, jpts 
    916                 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)  
     942            DO jn =1, jpts 
     943 
     944               IF ( Agrif_Parent(ln_zps) ) THEN ! Account for partial cells  
     945                                                ! linear vertical interpolation 
     946                  DO jj=j1,j2 
     947                     DO ji=i1,i2 
     948                        ! 
     949                        N_in  = mbkt(ji,jj) 
     950                        N_out = mbkt(ji,jj) 
     951                        z_in(1) = ptab(ji,jj,1,n2) 
     952                        tabin(1,1:jpts) = ptab(ji,jj,1,1:jpts) 
     953                        DO jk=2, N_in 
     954                           z_in(jk) = ptab(ji,jj,jk,n2) 
     955                           tabin(jk,1:jpts) = ptab(ji,jj,jk,1:jpts) 
     956                        END DO 
     957                        IF (.NOT.ln_linssh) z_in(1:N_in) = z_in(1:N_in) - ptab(ji,jj,k2,n2) 
     958                        z_out(1) = 0.5_wp * e3t(ji,jj,1,Krhs_a) 
     959                        DO jk=2, N_out 
     960                           z_out(jk) = z_out(jk-1) + 0.5_wp * (e3t(ji,jj,jk-1,Krhs_a) + e3t(ji,jj,jk,Krhs_a)) 
     961                        END DO 
     962                        IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out) - ssh(ji,jj,Krhs_a) 
     963                        CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),ptab(ji,jj,1:N_out,1:jpts), & 
     964                                      &   z_out(1:N_out),N_in,N_out,jpts)   
     965                     END DO 
     966                  END DO 
     967 
     968               ENDIF 
     969 
     970               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) 
     971 
    917972            END DO 
    918973         ENDIF 
     
    10201075                  zhtot = 0._wp 
    10211076                  DO jk=1,N_in 
    1022                      IF (jk==N_in) THEN 
    1023                         h_in(jk) = hu0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot 
    1024                      ELSE 
    1025                         h_in(jk) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy)  
    1026                      ENDIF 
     1077                     !IF (jk==N_in) THEN 
     1078                     !   h_in(jk) = hu0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot 
     1079                     !ELSE 
     1080                     !   h_in(jk) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy)  
     1081                     !ENDIF 
     1082                     h_in(jk) = e3u0_parent(ji,jj,jk)  
    10271083                     zhtot = zhtot + h_in(jk) 
    10281084                     IF( h_in(jk) .GT. 0. ) THEN 
     
    11421198                  zhtot = 0._wp 
    11431199                  DO jk=1,N_in 
    1144                      IF (jk==N_in) THEN 
    1145                         h_in(jk) = hv0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot 
    1146                      ELSE 
    1147                         h_in(jk) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox)  
    1148                      ENDIF 
     1200                     !IF (jk==N_in) THEN 
     1201                     !   h_in(jk) = hv0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot 
     1202                     !ELSE 
     1203                     !   h_in(jk) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox)  
     1204                     !ENDIF 
     1205                     h_in(jk) = e3v0_parent(ji,jj,jk) 
    11491206                     zhtot = zhtot + h_in(jk) 
    11501207                     IF( h_in(jk) .GT. 0. ) THEN 
     
    13261383   END SUBROUTINE interpub2b 
    13271384    
     1385   SUBROUTINE interpub2b_const( ptab, i1, i2, j1, j2, before ) 
     1386      !!---------------------------------------------------------------------- 
     1387      !!                  ***  ROUTINE interpub2b_const  *** 
     1388      !!----------------------------------------------------------------------   
     1389      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     1390      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     1391      LOGICAL                         , INTENT(in   ) ::   before 
     1392      ! 
     1393      REAL(wp) :: zrhoy 
     1394      !!----------------------------------------------------------------------   
     1395      IF( before ) THEN 
     1396         IF ( ln_bt_fw ) THEN 
     1397            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2) 
     1398         ELSE 
     1399            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * un_adv(i1:i2,j1:j2) 
     1400         ENDIF 
     1401      ELSE 
     1402         zrhoy = Agrif_Rhoy() 
     1403         ! 
     1404         ubdy(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) &  
     1405                           & / (zrhoy*e2u(i1:i2,j1:j2)) * umask(i1:i2,j1:j2,1) 
     1406         ! 
     1407      ENDIF 
     1408      !  
     1409   END SUBROUTINE interpub2b_const 
     1410 
     1411 
     1412   SUBROUTINE ub2b_cor( ptab, i1, i2, j1, j2, before ) 
     1413      !!---------------------------------------------------------------------- 
     1414      !!                  ***  ROUTINE ub2b_cor  *** 
     1415      !!----------------------------------------------------------------------   
     1416      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     1417      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     1418      LOGICAL                         , INTENT(in   ) ::   before 
     1419      ! 
     1420      INTEGER  :: ji, jj 
     1421      REAL(wp) :: zrhox, zrhoy, zx 
     1422      !!----------------------------------------------------------------------   
     1423      IF( before ) THEN 
     1424         ptab(:,:) = 0._wp 
     1425         DO ji=i1+1,i2-1 
     1426            DO jj=j1+1,j2 
     1427               ptab(ji,jj) = 0.25_wp*( ( vb2_b(ji+1,jj  )*e1v(ji+1,jj  )   &  
     1428                           &            -vb2_b(ji-1,jj  )*e1v(ji-1,jj  ) ) & 
     1429                           &          -( vb2_b(ji+1,jj-1)*e1v(ji+1,jj-1)   & 
     1430                           &            -vb2_b(ji-1,jj-1)*e1v(ji-1,jj-1) ) ) 
     1431            END DO 
     1432         END DO  
     1433      ELSE 
     1434         ! 
     1435         zrhox = Agrif_Rhox()  
     1436         zrhoy = Agrif_Rhoy() 
     1437         DO ji=i1,i2 
     1438            DO jj=j1,j2 
     1439               IF (utint_stage(ji,jj)==0) THEN  
     1440                  zx = 2._wp*MOD(ABS(mig0(ji)-nbghostcells-1), INT(Agrif_Rhox()))/zrhox - 1._wp   
     1441                  ubdy(ji,jj) = ubdy(ji,jj) + 0.25_wp*(1._wp-zx*zx) * ptab(ji,jj) &  
     1442                              &         / zrhoy *r1_e2u(ji,jj) * umask(ji,jj,1)  
     1443                  utint_stage(ji,jj) = 1  
     1444               ENDIF 
     1445            END DO 
     1446         END DO  
     1447         ! 
     1448      ENDIF 
     1449      !  
     1450   END SUBROUTINE ub2b_cor 
     1451 
    13281452 
    13291453   SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before ) 
     
    13631487 
    13641488 
     1489   SUBROUTINE interpvb2b_const( ptab, i1, i2, j1, j2, before ) 
     1490      !!---------------------------------------------------------------------- 
     1491      !!                  ***  ROUTINE interpub2b_const  *** 
     1492      !!----------------------------------------------------------------------   
     1493      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     1494      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     1495      LOGICAL                         , INTENT(in   ) ::   before 
     1496      ! 
     1497      REAL(wp) :: zrhox 
     1498      !!----------------------------------------------------------------------   
     1499      IF( before ) THEN 
     1500         IF ( ln_bt_fw ) THEN 
     1501            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2) 
     1502         ELSE 
     1503            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2) 
     1504         ENDIF 
     1505      ELSE 
     1506         zrhox = Agrif_Rhox() 
     1507         ! 
     1508         vbdy(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) & 
     1509                           & / (zrhox*e1v(i1:i2,j1:j2)) * vmask(i1:i2,j1:j2,1) 
     1510         ! 
     1511      ENDIF 
     1512      !  
     1513   END SUBROUTINE interpvb2b_const 
     1514 
     1515  
     1516   SUBROUTINE vb2b_cor( ptab, i1, i2, j1, j2, before ) 
     1517      !!---------------------------------------------------------------------- 
     1518      !!                  ***  ROUTINE vb2b_cor  *** 
     1519      !!----------------------------------------------------------------------   
     1520      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     1521      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     1522      LOGICAL                         , INTENT(in   ) ::   before 
     1523      ! 
     1524      INTEGER  :: ji, jj 
     1525      REAL(wp) :: zrhox, zrhoy, zy 
     1526      !!----------------------------------------------------------------------   
     1527      IF( before ) THEN 
     1528         ptab(:,:) = 0._wp 
     1529         DO ji=i1+1,i2-1 
     1530            DO jj=j1+1,j2 
     1531               ptab(ji,jj) = 0.25_wp*( ( ub2_b(ji  ,jj+1)*e2u(ji  ,jj+1)   &  
     1532                           &            -ub2_b(ji  ,jj-1)*e2u(ji  ,jj-1) ) & 
     1533                           &          -( ub2_b(ji-1,jj+1)*e2u(ji-1,jj+1)   & 
     1534                           &            -ub2_b(ji-1,jj-1)*e2u(ji-1,jj-1) ) ) 
     1535            END DO 
     1536         END DO  
     1537      ELSE 
     1538         ! 
     1539         zrhox = Agrif_Rhox()  
     1540         zrhoy = Agrif_Rhoy() 
     1541         DO ji=i1,i2 
     1542            DO jj=j1,j2 
     1543               IF (vtint_stage(ji,jj)==0) THEN  
     1544                  zy = 2._wp*MOD(ABS(mjg0(jj)-nbghostcells-1), INT(Agrif_Rhoy()))/zrhoy - 1._wp   
     1545                  vbdy(ji,jj) = vbdy(ji,jj) + 0.25_wp*(1._wp-zy*zy) * ptab(ji,jj) &  
     1546                              &         / zrhox * r1_e1v(ji,jj) * vmask(ji,jj,1)  
     1547                  vtint_stage(ji,jj) = 1  
     1548               ENDIF 
     1549            END DO 
     1550         END DO  
     1551         ! 
     1552      ENDIF 
     1553      !  
     1554   END SUBROUTINE vb2b_cor 
     1555 
     1556 
    13651557   SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before ) 
    13661558      !!---------------------------------------------------------------------- 
     
    13941586      !  
    13951587   END SUBROUTINE interpe3t 
     1588 
     1589 
     1590   SUBROUTINE interpe3t0_vremap( ptab, i1, i2, j1, j2, k1, k2, before ) 
     1591      !!---------------------------------------------------------------------- 
     1592      !!                  ***  ROUTINE interpe3t0_vremap  *** 
     1593      !!----------------------------------------------------------------------   
     1594      INTEGER                              , INTENT(in   ) :: i1, i2, j1, j2, k1, k2 
     1595      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab 
     1596      LOGICAL                              , INTENT(in   ) :: before 
     1597      ! 
     1598      INTEGER :: ji, jj, jk 
     1599      REAL(wp) :: zh 
     1600      !!----------------------------------------------------------------------   
     1601      !     
     1602      IF( before ) THEN 
     1603         IF ( ln_zps ) THEN 
     1604            DO jk = k1, k2 
     1605               DO jj = j1, j2 
     1606                  DO ji = i1, i2 
     1607                     ptab(ji, jj, jk) = e3t_1d(jk) 
     1608                  END DO 
     1609               END DO 
     1610            END DO 
     1611         ELSE 
     1612            ptab(i1:i2,j1:j2,k1:k2) = e3t_0(i1:i2,j1:j2,k1:k2) 
     1613         ENDIF 
     1614      ELSE 
     1615         ! 
     1616         DO jk = k1, k2 
     1617            DO jj = j1, j2 
     1618               DO ji = i1, i2 
     1619                  e3t0_parent(ji,jj,jk) = ptab(ji,jj,jk) 
     1620               END DO 
     1621            END DO 
     1622         END DO 
     1623 
     1624         ! Retrieve correct scale factor at the bottom: 
     1625         DO jj = j1, j2 
     1626            DO ji = i1, i2 
     1627               zh = 0._wp 
     1628               DO jk = 1, mbkt_parent(ji, jj)-1 
     1629                  zh = zh + e3t0_parent(ji,jj,jk) 
     1630               END DO 
     1631               e3t0_parent(ji,jj,mbkt_parent(ji,jj)) = ht0_parent(ji, jj) - zh 
     1632            END DO 
     1633         END DO 
     1634          
     1635      ENDIF 
     1636      !  
     1637   END SUBROUTINE interpe3t0_vremap 
     1638 
    13961639 
    13971640   SUBROUTINE interpglamt( ptab, i1, i2, j1, j2, before ) 
  • NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_oce_sponge.F90

    r13565 r13937  
    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 
     
    4546      !!---------------------------------------------------------------------- 
    4647      REAL(wp) ::   zcoef   ! local scalar 
    47        
     48      INTEGER  :: istart, iend, jstart, jend  
    4849      !!---------------------------------------------------------------------- 
    4950      ! 
     
    5253      zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot()) 
    5354 
    54       CALL Agrif_Sponge 
    5555      Agrif_SpecialValue    = 0._wp 
    5656      Agrif_UseSpecialValue = .TRUE. 
     
    9393      tabspongedone_v = .FALSE. 
    9494      CALL Agrif_Bc_Variable( vn_sponge_id, calledweight=zcoef, procname=interpvn_sponge ) 
     95       
     96      IF ( nn_shift_bar>0 ) THEN ! then split sponge between 2d and 3d 
     97         zcoef = REAL(Agrif_NbStepint(),wp)/REAL(Agrif_rhot()) ! forward tsplit 
     98         tabspongedone_u = .FALSE. 
     99         tabspongedone_v = .FALSE.          
     100         CALL Agrif_Bc_Variable( unb_sponge_id, calledweight=zcoef, procname=interpunb_sponge ) 
     101         ! 
     102         tabspongedone_u = .FALSE. 
     103         tabspongedone_v = .FALSE. 
     104         CALL Agrif_Bc_Variable( vnb_sponge_id, calledweight=zcoef, procname=interpvnb_sponge ) 
     105      ENDIF 
    95106      ! 
    96107      Agrif_UseSpecialValue = .FALSE. 
    97108      use_sign_north        = .FALSE. 
    98109      l_vremap              = .FALSE. 
     110      ! 
    99111#endif 
    100112      ! 
     
    127139 
    128140#if defined SPONGE || defined SPONGE_TOP 
    129       IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN 
    130          ! Define ramp from boundaries towards domain interior at F-points 
    131          ! Store it in ztabramp 
    132  
    133          ispongearea    = nn_sponge_len * Agrif_irhox() 
    134          z1_ispongearea = 1._wp / REAL( MAX(ispongearea,1), wp ) 
    135          jspongearea    = nn_sponge_len * Agrif_irhoy() 
    136          z1_jspongearea = 1._wp / REAL( MAX(jspongearea,1), wp ) 
     141      ! Define ramp from boundaries towards domain interior at F-points 
     142      ! Store it in ztabramp 
     143 
     144      ispongearea    = nn_sponge_len * Agrif_irhox() 
     145      z1_ispongearea = 1._wp / REAL( MAX(ispongearea,1), wp ) 
     146      jspongearea    = nn_sponge_len * Agrif_irhoy() 
     147      z1_jspongearea = 1._wp / REAL( MAX(jspongearea,1), wp ) 
    137148          
    138          ztabramp(:,:) = 0._wp 
    139  
    140          IF( lk_west ) THEN                             ! --- West --- ! 
    141             ind1 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
    142             ind2 = nn_hls + 1 + nbghostcells + ispongearea  
    143             DO ji = mi0(ind1), mi1(ind2)    
    144                DO jj = 1, jpj                
    145                   ztabramp(ji,jj) =                       REAL(ind2 - mig(ji), wp) * z1_ispongearea 
    146                END DO 
    147             END DO 
    148             ! ghost cells: 
    149             ind1 = 1 
    150             ind2 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
    151             DO ji = mi0(ind1), mi1(ind2)    
    152                DO jj = 1, jpj                
    153                   ztabramp(ji,jj) = 1._wp 
    154                END DO 
    155             END DO 
    156          ENDIF 
    157          IF( lk_east ) THEN                             ! --- East --- ! 
    158             ind1 = jpiglo - ( nn_hls + nbghostcells ) - ispongearea - 1 
    159             ind2 = jpiglo - ( nn_hls + nbghostcells ) - 1    ! halo + land + nbghostcells - 1 
    160             DO ji = mi0(ind1), mi1(ind2) 
    161                DO jj = 1, jpj 
    162                   ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mig(ji) - ind1, wp) * z1_ispongearea )  
    163                END DO 
    164             END DO 
    165             ! ghost cells: 
    166             ind1 = jpiglo - ( nn_hls + nbghostcells ) - 1    ! halo + land + nbghostcells - 1 
    167             ind2 = jpiglo - 1 
    168             DO ji = mi0(ind1), mi1(ind2) 
    169                DO jj = 1, jpj 
    170                   ztabramp(ji,jj) = 1._wp 
    171                END DO 
    172             END DO 
    173          ENDIF       
    174          IF( lk_south ) THEN                            ! --- South --- ! 
    175             ind1 = nn_hls + 1 + nbghostcells                 ! halo + land + nbghostcells 
    176             ind2 = nn_hls + 1 + nbghostcells + jspongearea  
    177             DO jj = mj0(ind1), mj1(ind2)  
    178                DO ji = 1, jpi 
    179                   ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(ind2 - mjg(jj), wp) * z1_jspongearea ) 
    180                END DO 
    181             END DO 
    182             ! ghost cells: 
    183             ind1 = 1 
    184             ind2 = nn_hls + 1 + nbghostcells                 ! halo + land + nbghostcells 
    185             DO jj = mj0(ind1), mj1(ind2)  
    186                DO ji = 1, jpi 
    187                   ztabramp(ji,jj) = 1._wp 
    188                END DO 
    189             END DO 
    190          ENDIF 
    191          IF( lk_north ) THEN                            ! --- North --- ! 
    192             ind1 = jpjglo - ( nn_hls + nbghostcells ) - jspongearea - 1 
    193             ind2 = jpjglo - ( nn_hls + nbghostcells ) - 1    ! halo + land + nbghostcells - 1 
    194             DO jj = mj0(ind1), mj1(ind2) 
    195                DO ji = 1, jpi 
    196                   ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mjg(jj) - ind1, wp) * z1_jspongearea )  
    197                END DO 
    198             END DO 
    199             ! ghost cells: 
    200             ind1 = jpjglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    201             ind2 = jpjglo 
    202             DO jj = mj0(ind1), mj1(ind2) 
    203                DO ji = 1, jpi 
    204                   ztabramp(ji,jj) = 1._wp 
    205                END DO 
    206             END DO 
    207          ENDIF 
    208          ! 
     149      ztabramp(:,:) = 0._wp 
     150 
     151      IF( lk_west ) THEN                             ! --- West --- ! 
     152         ind1 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
     153         ind2 = nn_hls + 1 + nbghostcells + ispongearea  
     154         DO ji = mi0(ind1), mi1(ind2)    
     155            DO jj = 1, jpj                
     156               ztabramp(ji,jj) =                       REAL(ind2 - mig(ji), wp) * z1_ispongearea 
     157            END DO 
     158         END DO 
     159         ! ghost cells: 
     160         ind1 = 1 
     161         ind2 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
     162         DO ji = mi0(ind1), mi1(ind2)    
     163            DO jj = 1, jpj                
     164               ztabramp(ji,jj) = 1._wp 
     165            END DO 
     166         END DO 
    209167      ENDIF 
    210  
     168      IF( lk_east ) THEN                             ! --- East --- ! 
     169         ind1 = jpiglo - ( nn_hls + nbghostcells ) - ispongearea - 1 
     170         ind2 = jpiglo - ( nn_hls + nbghostcells ) - 1    ! halo + land + nbghostcells - 1 
     171         DO ji = mi0(ind1), mi1(ind2) 
     172            DO jj = 1, jpj 
     173               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mig(ji) - ind1, wp) * z1_ispongearea )  
     174            END DO 
     175         END DO 
     176         ! ghost cells: 
     177         ind1 = jpiglo - ( nn_hls + nbghostcells ) - 1    ! halo + land + nbghostcells - 1 
     178         ind2 = jpiglo - 1 
     179         DO ji = mi0(ind1), mi1(ind2) 
     180            DO jj = 1, jpj 
     181               ztabramp(ji,jj) = 1._wp 
     182            END DO 
     183         END DO 
     184      ENDIF       
     185      IF( lk_south ) THEN                            ! --- South --- ! 
     186         ind1 = nn_hls + 1 + nbghostcells                 ! halo + land + nbghostcells 
     187         ind2 = nn_hls + 1 + nbghostcells + jspongearea  
     188         DO jj = mj0(ind1), mj1(ind2)  
     189            DO ji = 1, jpi 
     190               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(ind2 - mjg(jj), wp) * z1_jspongearea ) 
     191            END DO 
     192         END DO 
     193         ! ghost cells: 
     194         ind1 = 1 
     195         ind2 = nn_hls + 1 + nbghostcells                 ! halo + land + nbghostcells 
     196         DO jj = mj0(ind1), mj1(ind2)  
     197            DO ji = 1, jpi 
     198               ztabramp(ji,jj) = 1._wp 
     199            END DO 
     200         END DO 
     201      ENDIF 
     202      IF( lk_north ) THEN                            ! --- North --- ! 
     203         ind1 = jpjglo - ( nn_hls + nbghostcells ) - jspongearea - 1 
     204         ind2 = jpjglo - ( nn_hls + nbghostcells ) - 1    ! halo + land + nbghostcells - 1 
     205         DO jj = mj0(ind1), mj1(ind2) 
     206            DO ji = 1, jpi 
     207               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mjg(jj) - ind1, wp) * z1_jspongearea )  
     208            END DO 
     209         END DO 
     210         ! ghost cells: 
     211         ind1 = jpjglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
     212         ind2 = jpjglo 
     213         DO jj = mj0(ind1), mj1(ind2) 
     214            DO ji = 1, jpi 
     215               ztabramp(ji,jj) = 1._wp 
     216            END DO 
     217         END DO 
     218      ENDIF       
     219      ! 
    211220      ! Tracers 
    212       IF( .NOT. spongedoneT ) THEN 
    213          fspu(:,:) = 0._wp 
    214          fspv(:,:) = 0._wp 
    215          DO_2D( 0, 0, 0, 0 ) 
    216             fspu(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji,jj-1) ) * ssumask(ji,jj) 
    217             fspv(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji-1,jj) ) * ssvmask(ji,jj) 
    218          END_2D 
    219       ENDIF 
     221      fspu(:,:) = 0._wp 
     222      fspv(:,:) = 0._wp 
     223      DO_2D( 0, 0, 0, 0 ) 
     224         fspu(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji,jj-1) ) * ssumask(ji,jj) 
     225         fspv(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji-1,jj) ) * ssvmask(ji,jj) 
     226      END_2D 
    220227 
    221228      ! Dynamics 
    222       IF( .NOT. spongedoneU ) THEN 
    223          fspt(:,:) = 0._wp 
    224          fspf(:,:) = 0._wp 
    225          DO_2D( 0, 0, 0, 0 ) 
    226             fspt(ji,jj) = 0.25_wp * ( ztabramp(ji  ,jj  ) + ztabramp(ji-1,jj  ) & 
    227                                   &  +ztabramp(ji  ,jj-1) + ztabramp(ji-1,jj-1) ) * ssmask(ji,jj) 
    228             fspf(ji,jj) = ztabramp(ji,jj) * ssvmask(ji,jj) * ssvmask(ji,jj+1) 
    229          END_2D 
    230       ENDIF 
     229      fspt(:,:) = 0._wp 
     230      fspf(:,:) = 0._wp 
     231      DO_2D( 0, 0, 0, 0 ) 
     232         fspt(ji,jj) = 0.25_wp * ( ztabramp(ji  ,jj  ) + ztabramp(ji-1,jj  ) & 
     233                               &  +ztabramp(ji  ,jj-1) + ztabramp(ji-1,jj-1) ) * ssmask(ji,jj) 
     234         fspf(ji,jj) = ztabramp(ji,jj) * ssvmask(ji,jj) * ssvmask(ji,jj+1) 
     235      END_2D 
    231236       
    232       IF( .NOT. spongedoneT .AND. .NOT. spongedoneU ) THEN 
    233          CALL lbc_lnk_multi( 'agrif_Sponge', fspu, 'U', 1._wp, fspv, 'V', 1._wp, fspt, 'T', 1._wp, fspf, 'F', 1._wp ) 
    234          spongedoneT = .TRUE. 
    235          spongedoneU = .TRUE. 
    236       ENDIF 
    237       IF( .NOT. spongedoneT ) THEN 
    238          CALL lbc_lnk_multi( 'agrif_Sponge', fspu, 'U', 1._wp, fspv, 'V', 1._wp ) 
    239          spongedoneT = .TRUE. 
    240       ENDIF 
    241       IF( .NOT. spongedoneT .AND. .NOT. spongedoneU ) THEN 
    242          CALL lbc_lnk_multi( 'agrif_Sponge', fspt, 'T', 1._wp, fspf, 'F', 1._wp ) 
    243          spongedoneU = .TRUE. 
    244       ENDIF 
     237      CALL lbc_lnk_multi( 'agrif_Sponge', fspu, 'U', 1._wp, fspv, 'V', 1._wp, fspt, 'T', 1._wp, fspf, 'F', 1._wp ) 
    245238      ! 
    246239      ! Remove vertical interpolation where not needed: 
     
    254247   END SUBROUTINE Agrif_Sponge 
    255248 
    256     
    257    SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
     249 
     250   SUBROUTINE Agrif_Sponge_2d 
     251      !!---------------------------------------------------------------------- 
     252      !!                 *** ROUTINE  Agrif_Sponge_2d *** 
     253      !!---------------------------------------------------------------------- 
     254      INTEGER  ::   ji, jj, ind1, ind2, ishift, jshift 
     255      INTEGER  ::   ispongearea, jspongearea 
     256      REAL(wp) ::   z1_ispongearea, z1_jspongearea 
     257      REAL(wp), DIMENSION(jpi,jpj) :: ztabramp 
     258      !!---------------------------------------------------------------------- 
     259      ! 
     260      ! Sponge 1d example with: 
     261      !      iraf = 3 ; nbghost = 3 ; nn_sponge_len = 2 
     262      !                         
     263      !coarse :     U     T     U     T     U     T     U 
     264      !|            |           |           |           | 
     265      !fine :     t u t u t u t u t u t u t u t u t u t u t 
     266      !sponge val:0 1   1   1   1  5/6 4/6 3/6 2/6 1/6  0  
     267      !           |   ghost     | <-- sponge area  -- > | 
     268      !           |   points    |                       | 
     269      !                         |--> dynamical interface 
     270 
     271#if defined SPONGE || defined SPONGE_TOP 
     272      ! Define ramp from boundaries towards domain interior at F-points 
     273      ! Store it in ztabramp 
     274 
     275      ispongearea    = nn_sponge_len * Agrif_irhox() 
     276      z1_ispongearea = 1._wp / REAL( MAX(ispongearea,1), wp ) 
     277      jspongearea    = nn_sponge_len * Agrif_irhoy() 
     278      z1_jspongearea = 1._wp / REAL( MAX(jspongearea,1), wp ) 
     279      ishift = nn_shift_bar * Agrif_irhox() 
     280      jshift = nn_shift_bar * Agrif_irhoy() 
     281         
     282      ztabramp(:,:) = 0._wp 
     283 
     284      IF( lk_west ) THEN                             ! --- West --- ! 
     285         ind1 = nn_hls + 1 + nbghostcells + ishift 
     286         ind2 = nn_hls + 1 + nbghostcells + ishift + ispongearea  
     287         DO ji = mi0(ind1), mi1(ind2)    
     288            DO jj = 1, jpj                
     289               ztabramp(ji,jj) =                       REAL(ind2 - mig(ji), wp) * z1_ispongearea 
     290            END DO 
     291         END DO 
     292         ! ghost cells: 
     293         ind1 = 1 
     294         ind2 = nn_hls + 1 + nbghostcells + ishift               ! halo + land + nbghostcells 
     295         DO ji = mi0(ind1), mi1(ind2)    
     296            DO jj = 1, jpj                
     297               ztabramp(ji,jj) = 1._wp 
     298            END DO 
     299         END DO 
     300      ENDIF 
     301      IF( lk_east ) THEN                             ! --- East --- ! 
     302         ind1 = jpiglo - ( nn_hls + nbghostcells + ishift) - ispongearea - 1 
     303         ind2 = jpiglo - ( nn_hls + nbghostcells + ishift) - 1    ! halo + land + nbghostcells - 1 
     304         DO ji = mi0(ind1), mi1(ind2) 
     305            DO jj = 1, jpj 
     306               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mig(ji) - ind1, wp) * z1_ispongearea )  
     307            END DO 
     308         END DO 
     309         ! ghost cells: 
     310         ind1 = jpiglo - ( nn_hls + nbghostcells + ishift) - 1    ! halo + land + nbghostcells - 1 
     311         ind2 = jpiglo - 1 
     312         DO ji = mi0(ind1), mi1(ind2) 
     313            DO jj = 1, jpj 
     314               ztabramp(ji,jj) = 1._wp 
     315            END DO 
     316         END DO 
     317      ENDIF       
     318      IF( lk_south ) THEN                            ! --- South --- ! 
     319         ind1 = nn_hls + 1 + nbghostcells + jshift                ! halo + land + nbghostcells 
     320         ind2 = nn_hls + 1 + nbghostcells + jshift + jspongearea  
     321         DO jj = mj0(ind1), mj1(ind2)  
     322            DO ji = 1, jpi 
     323               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(ind2 - mjg(jj), wp) * z1_jspongearea ) 
     324            END DO 
     325         END DO 
     326         ! ghost cells: 
     327         ind1 = 1 
     328         ind2 = nn_hls + 1 + nbghostcells + jshift                ! halo + land + nbghostcells 
     329         DO jj = mj0(ind1), mj1(ind2)  
     330            DO ji = 1, jpi 
     331               ztabramp(ji,jj) = 1._wp 
     332            END DO 
     333         END DO 
     334      ENDIF 
     335      IF( lk_north ) THEN                            ! --- North --- ! 
     336         ind1 = jpjglo - ( nn_hls + nbghostcells + jshift) - jspongearea - 1 
     337         ind2 = jpjglo - ( nn_hls + nbghostcells + jshift) - 1    ! halo + land + nbghostcells - 1 
     338         DO jj = mj0(ind1), mj1(ind2) 
     339            DO ji = 1, jpi 
     340               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mjg(jj) - ind1, wp) * z1_jspongearea )  
     341            END DO 
     342         END DO 
     343         ! ghost cells: 
     344         ind1 = jpjglo - ( nn_hls + nbghostcells + jshift)      ! halo + land + nbghostcells - 1 
     345         ind2 = jpjglo 
     346         DO jj = mj0(ind1), mj1(ind2) 
     347            DO ji = 1, jpi 
     348               ztabramp(ji,jj) = 1._wp 
     349            END DO 
     350         END DO 
     351      ENDIF 
     352      ! 
     353      ! Tracers 
     354      fspu_2d(:,:) = 0._wp 
     355      fspv_2d(:,:) = 0._wp 
     356      DO_2D( 0, 0, 0, 0 ) 
     357         fspu_2d(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji,jj-1) ) * ssumask(ji,jj) 
     358         fspv_2d(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji-1,jj) ) * ssvmask(ji,jj) 
     359      END_2D 
     360 
     361      ! Dynamics 
     362      fspt_2d(:,:) = 0._wp 
     363      fspf_2d(:,:) = 0._wp 
     364      DO_2D( 0, 0, 0, 0 ) 
     365         fspt_2d(ji,jj) = 0.25_wp * ( ztabramp(ji  ,jj  ) + ztabramp(ji-1,jj  ) & 
     366                               &  +ztabramp(ji  ,jj-1) + ztabramp(ji-1,jj-1) ) * ssmask(ji,jj) 
     367         fspf_2d(ji,jj) = ztabramp(ji,jj) * ssvmask(ji,jj) * ssvmask(ji,jj+1) 
     368         END_2D 
     369      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 ) 
     370      ! 
     371#endif 
     372      ! 
     373   END SUBROUTINE Agrif_Sponge_2d 
     374 
     375 
     376   SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before)  
    258377      !!---------------------------------------------------------------------- 
    259378      !!                 *** ROUTINE interptsn_sponge *** 
     
    266385      INTEGER  ::   iku, ikv 
    267386      REAL(wp) :: ztsa, zabe1, zabe2, zbtr, zhtot 
    268       REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv 
     387      REAL(wp), DIMENSION(i1-1:i2,j1-1:j2,jpk) :: ztu, ztv 
    269388      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff 
    270389      ! vertical interpolation: 
    271390      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child 
    272       REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin 
    273       REAL(wp), DIMENSION(k1:k2) :: h_in 
    274       REAL(wp), DIMENSION(1:jpk) :: h_out 
     391      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin, tabin_i 
     392      REAL(wp), DIMENSION(k1:k2) :: z_in, z_in_i, h_in_i 
     393      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out 
    275394      INTEGER :: N_in, N_out 
    276395      !!---------------------------------------------------------------------- 
     
    287406         END DO 
    288407 
    289         IF ( l_vremap ) THEN 
    290  
    291            ! Interpolate thicknesses 
    292            ! Warning: these are masked, hence extrapolated prior interpolation. 
    293            DO jk=k1,k2 
    294               DO jj=j1,j2 
    295                  DO ji=i1,i2 
    296                    tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kbb_a) 
    297                  END DO 
    298               END DO 
    299            END DO 
    300  
    301            ! Extrapolate thicknesses in partial bottom cells: 
    302            ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
    303            IF (ln_zps) THEN 
    304               DO jj=j1,j2 
    305                  DO ji=i1,i2 
    306                     jk = mbkt(ji,jj) 
    307                     tabres(ji,jj,jk,jpts+1) = 0._wp 
    308                  END DO 
    309               END DO            
    310            END IF 
    311       
    312            ! Save ssh at last level: 
    313            IF (.NOT.ln_linssh) THEN 
    314               tabres(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kbb_a)*tmask(i1:i2,j1:j2,1)  
    315            ELSE 
    316               tabres(i1:i2,j1:j2,k2,jpts+1) = 0._wp 
    317            END IF       
    318         END IF 
    319  
    320       ELSE    
    321          ! 
    322          IF ( l_vremap ) THEN 
    323  
    324             IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp 
    325  
     408         IF ( l_vremap.OR.ln_zps ) THEN 
     409 
     410            ! Fill cell depths (i.e. gdept) to be interpolated 
     411            ! Warning: these are masked, hence extrapolated prior interpolation. 
    326412            DO jj=j1,j2 
    327413               DO ji=i1,i2 
     414                  tabres(ji,jj,k1,jpts+1) = 0.5_wp * tmask(ji,jj,k1) * e3t(ji,jj,k1,Kbb_a) 
     415                  DO jk=k1+1,k2 
     416                     tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * & 
     417                        & ( tabres(ji,jj,jk-1,jpts+1) + 0.5_wp * (e3t(ji,jj,jk-1,Kbb_a)+e3t(ji,jj,jk,Kbb_a)) ) 
     418                  END DO 
     419               END DO 
     420            END DO 
     421 
     422            ! Save ssh at last level: 
     423            IF ( .NOT.ln_linssh ) THEN 
     424               tabres(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kbb_a)*tmask(i1:i2,j1:j2,1)  
     425            END IF   
     426     
     427         END IF 
     428 
     429      ELSE  
     430         ! 
     431         IF ( l_vremap ) THEN 
     432 
     433            IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp 
     434 
     435            DO jj=j1,j2 
     436               DO ji=i1,i2 
     437 
    328438                  tabres_child(ji,jj,:,:) = 0._wp  
     439                  ! Build vertical grids: 
    329440                  N_in = mbkt_parent(ji,jj) 
    330                   zhtot = 0._wp 
    331                   DO jk=1,N_in !k2 = jpk of parent grid 
    332                      IF (jk==N_in) THEN 
    333                         h_in(jk) = ht0_parent(ji,jj) + tabres(ji,jj,k2,n2) - zhtot 
    334                      ELSE 
    335                         h_in(jk) = tabres(ji,jj,jk,n2) 
    336                      END IF 
    337                      zhtot = zhtot + h_in(jk) 
    338                      tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1) 
    339                   END DO 
    340                   N_out = 0 
    341                   DO jk=1,jpk ! jpk of child grid 
    342                      IF (tmask(ji,jj,jk) == 0) EXIT  
    343                      N_out = N_out + 1 
    344                      h_out(jk) = e3t(ji,jj,jk,Kbb_a) !Child grid scale factors. Could multiply by e1e2t here instead of division above 
    345                   END DO 
     441                  ! Input grid (account for partial cells if any): 
     442                  DO jk=1,N_in 
     443                     z_in(jk) = tabres(ji,jj,jk,n2) - tabres(ji,jj,k2,n2) 
     444                     tabin(jk,1:jpts) = tabres(ji,jj,jk,1:jpts) 
     445                  END DO 
     446                   
     447                  ! Intermediate grid: 
     448                  DO jk = 1, N_in 
     449                     h_in_i(jk) = e3t0_parent(ji,jj,jk) * &  
     450                       &       (1._wp + tabres(ji,jj,k2,n2)/(ht0_parent(ji,jj)*ssmask(ji,jj) + 1._wp - ssmask(ji,jj))) 
     451                  END DO 
     452                  z_in_i(1) = 0.5_wp * h_in_i(1) 
     453                  DO jk=2,N_in 
     454                     z_in_i(jk) = z_in_i(jk-1) + 0.5_wp * ( h_in_i(jk) + h_in_i(jk-1) ) 
     455                  END DO 
     456                  z_in_i(1:N_in) = z_in_i(1:N_in)  - tabres(ji,jj,k2,n2) 
     457 
     458                  ! Output (Child) grid: 
     459                  N_out = mbkt(ji,jj) 
     460                  DO jk=1,N_out 
     461                     h_out(jk) = e3t(ji,jj,jk,Kbb_a) 
     462                  END DO 
     463                  z_out(1) = 0.5_wp * h_out(1) 
     464                  DO jk=2,N_out 
     465                     z_out(jk) = z_out(jk-1) + 0.5_wp * ( h_out(jk)+h_out(jk-1) ) 
     466                  END DO 
     467                  IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out)  - ssh(ji,jj,Kbb_a) 
    346468 
    347469                  ! Account for small differences in the free-surface 
    348                   IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 
    349                      h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 
     470                  IF ( sum(h_out(1:N_out)) > sum(h_in_i(1:N_in) )) THEN 
     471                     h_out(1) = h_out(1)  - ( sum(h_out(1:N_out))-sum(h_in_i(1:N_in)) ) 
    350472                  ELSE 
    351                      h_in(1)   = h_in(1)   - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 
     473                     h_in_i(1)= h_in_i(1) - ( sum(h_in_i(1:N_in))-sum(h_out(1:N_out)) ) 
    352474                  END IF 
    353475                  IF (N_in*N_out > 0) THEN 
    354                      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) 
     476                     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) 
     477                     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) 
     478!                     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)   
    355479                  ENDIF 
    356480               END DO 
     
    367491         ELSE 
    368492 
     493            IF ( Agrif_Parent(ln_zps) ) THEN ! Account for partial cells 
     494 
     495               DO jj=j1,j2 
     496                  DO ji=i1,i2 
     497                     ! 
     498                     N_in  = mbkt(ji,jj)  
     499                     N_out = mbkt(ji,jj)  
     500                     z_in(1) = tabres(ji,jj,1,n2) 
     501                     tabin(1,1:jpts) = tabres(ji,jj,1,1:jpts) 
     502                     DO jk=2, N_in 
     503                        z_in(jk) = tabres(ji,jj,jk,n2) 
     504                        tabin(jk,1:jpts) = tabres(ji,jj,jk,1:jpts) 
     505                     END DO  
     506                     IF (.NOT.ln_linssh) z_in(1:N_in) = z_in(1:N_in) - tabres(ji,jj,k2,n2) 
     507 
     508                     z_out(1) = 0.5_wp * e3t(ji,jj,1,Kbb_a) 
     509                     DO jk=2, N_out 
     510                        z_out(jk) = z_out(jk-1) + 0.5_wp * (e3t(ji,jj,jk-1,Kbb_a) + e3t(ji,jj,jk,Kbb_a))  
     511                     END DO  
     512                     IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out) - ssh(ji,jj,Kbb_a) 
     513 
     514                     CALL remap_linear(tabin(1:N_in,1:jpts), z_in(1:N_in), tabres(ji,jj,1:N_out,1:jpts), & 
     515                                         &   z_out(1:N_out), N_in, N_out, jpts) 
     516                  END DO 
     517               END DO 
     518            ENDIF 
     519 
    369520            DO jj=j1,j2 
    370521               DO ji=i1,i2 
     
    374525               END DO 
    375526            END DO 
     527 
    376528         END IF 
    377529 
    378530         DO jn = 1, jpts             
    379531            DO jk = 1, jpkm1 
    380                ztu(i1:i2,j1:j2,jk) = 0._wp 
     532               ztu(i1-1:i2,j1-1:j2,jk) = 0._wp 
    381533               DO jj = j1,j2 
    382534                  DO ji = i1,i2-1 
    383                      zabe1 = fspu(ji,jj) * rn_sponge_tra * r1_Dt * umask(ji,jj,jk) * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) 
    384                      ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) )  
    385                   END DO 
    386                END DO 
    387                ztv(i1:i2,j1:j2,jk) = 0._wp 
     535                     zabe1 = rn_sponge_tra * r1_Dt * umask(ji,jj,jk) * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) 
     536                     ztu(ji,jj,jk) = zabe1 * fspu(ji,jj) * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) )  
     537                  END DO 
     538               END DO 
     539               ztv(i1-1:i2,j1-1:j2,jk) = 0._wp 
    388540               DO ji = i1,i2 
    389541                  DO jj = j1,j2-1 
    390                      zabe2 = fspv(ji,jj) * rn_sponge_tra * r1_Dt * vmask(ji,jj,jk) * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a) 
    391                      ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
     542                     zabe2 = rn_sponge_tra * r1_Dt * vmask(ji,jj,jk) * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a) 
     543                     ztv(ji,jj,jk) = zabe2 * fspv(ji,jj) * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
    392544                  END DO 
    393545               END DO 
     
    406558            END DO 
    407559            ! 
     560! JC: there is something wrong with the Laplacian in corners 
    408561            DO jk = 1, jpkm1 
    409                DO jj = j1+1,j2-1 
    410                   DO ji = i1+1,i2-1 
     562               DO jj = j1,j2 
     563                  DO ji = i1,i2 
    411564                     IF (.NOT. tabspongedone_tsn(ji,jj)) THEN  
    412565                        zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a) 
    413566                        ! horizontal diffusive trends 
    414                         ztsa = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) & 
    415                              &  - rn_trelax_tra * r1_Dt * fspt(ji,jj) * tsbdiff(ji,jj,jk,jn)  
     567                        ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   &  
     568                          &           + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) & 
     569                          &   - rn_trelax_tra * r1_Dt * fspt(ji,jj) * tsbdiff(ji,jj,jk,jn) 
    416570                        ! add it to the general tracer trends 
    417571                        ts(ji,jj,jk,jn,Krhs_a) = ts(ji,jj,jk,jn,Krhs_a) + ztsa 
     
    419573                  END DO 
    420574               END DO 
     575 
    421576            END DO 
    422577            ! 
    423578         END DO 
    424579         ! 
    425          tabspongedone_tsn(i1+1:i2-1,j1+1:j2-1) = .TRUE. 
     580         tabspongedone_tsn(i1:i2,j1:j2) = .TRUE. 
    426581         ! 
    427582      ENDIF 
     
    503658                  zhtot = 0._wp 
    504659                  DO jk=1,N_in 
    505                      IF (jk==N_in) THEN 
    506                         h_in(jk) = hu0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot 
    507                      ELSE 
    508                         h_in(jk) = tabres(ji,jj,jk,m2) 
    509                      ENDIF 
     660                     !IF (jk==N_in) THEN 
     661                     !   h_in(jk) = hu0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot 
     662                     !ELSE 
     663                     !   h_in(jk) = tabres(ji,jj,jk,m2) 
     664                     !ENDIF 
     665                     h_in(jk) = e3u0_parent(ji,jj,jk) 
    510666                     zhtot = zhtot + h_in(jk) 
    511667                     tabin(jk) = tabres(ji,jj,jk,m1) 
     
    533689 
    534690            ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:) 
    535  
    536691         ELSE 
    537692 
    538693            ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:) 
    539  
     694   
    540695         ENDIF 
    541696         ! 
     
    688843                  zhtot = 0._wp 
    689844                  DO jk=1,N_in 
    690                      IF (jk==N_in) THEN 
    691                         h_in(jk) = hv0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot 
    692                      ELSE 
    693                         h_in(jk) = tabres(ji,jj,jk,m2) 
    694                      ENDIF 
     845                     !IF (jk==N_in) THEN 
     846                     !   h_in(jk) = hv0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot 
     847                     !ELSE 
     848                     !   h_in(jk) = tabres(ji,jj,jk,m2) 
     849                     !ENDIF 
     850                     h_in(jk) = e3v0_parent(ji,jj,jk) 
    695851                     zhtot = zhtot + h_in(jk) 
    696852                     tabin(jk) = tabres(ji,jj,jk,m1) 
     
    718874 
    719875            vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:)   
    720  
    721876         ELSE 
    722877 
     
    786941      ! 
    787942   END SUBROUTINE interpvn_sponge 
     943 
     944   SUBROUTINE interpunb_sponge(tabres,i1,i2,j1,j2, before) 
     945      !!--------------------------------------------- 
     946      !!   *** ROUTINE interpunb_sponge *** 
     947      !!---------------------------------------------     
     948      INTEGER, INTENT(in) :: i1,i2,j1,j2 
     949      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres 
     950      LOGICAL, INTENT(in) :: before 
     951 
     952      INTEGER  :: ji, jj, ind1, jmax 
     953      ! sponge parameters  
     954      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr 
     955      REAL(wp), DIMENSION(i1:i2,j1:j2) :: ubdiff 
     956      REAL(wp), DIMENSION(i1:i2,j1:j2) :: rotdiff, hdivdiff 
     957      !!---------------------------------------------     
     958      ! 
     959      IF( before ) THEN 
     960         DO jj=j1,j2 
     961            DO ji=i1,i2 
     962               tabres(ji,jj) = uu_b(ji,jj,Kmm_a) 
     963            END DO 
     964         END DO 
     965 
     966      ELSE 
     967 
     968         ubdiff(i1:i2,j1:j2) = (uu_b(i1:i2,j1:j2,Kmm_a) - tabres(i1:i2,j1:j2))*umask(i1:i2,j1:j2,1) 
     969         ! 
     970         !                                             ! -------- 
     971         ! Horizontal divergence                       !   div 
     972         !                                             ! -------- 
     973         DO jj = j1,j2 
     974            DO ji = i1+1,i2   ! vector opt. 
     975               zbtr = rn_sponge_dyn * r1_Dt * fspt_2d(ji,jj) * r1_ht_0(ji,jj) 
     976               hdivdiff(ji,jj) = (  e2u(ji  ,jj)*hu(ji  ,jj,Kbb_a) * ubdiff(ji  ,jj) & 
     977                                  &-e2u(ji-1,jj)*hu(ji-1,jj,Kbb_a) * ubdiff(ji-1,jj) ) * zbtr 
     978            END DO 
     979         END DO 
     980 
     981         DO jj = j1,j2-1 
     982            DO ji = i1,i2   ! vector opt. 
     983               zbtr = rn_sponge_dyn * r1_Dt * fspf_2d(ji,jj) * hf_0(ji,jj) 
     984               rotdiff(ji,jj) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1)   & 
     985                              &   +e1u(ji,jj  ) * ubdiff(ji,jj  ) ) * fmask(ji,jj,1) * zbtr  
     986            END DO 
     987         END DO 
     988         ! 
     989         DO jj = j1+1, j2-1 
     990            DO ji = i1+1, i2-1   ! vector opt. 
     991               IF (.NOT. tabspongedone_u(ji,jj)) THEN 
     992                  ze2u = rotdiff (ji,jj) 
     993                  ze1v = hdivdiff(ji,jj) 
     994                  ! horizontal diffusive trends 
     995                  zua = - ( ze2u - rotdiff (ji,jj-1) ) * r1_e2u(ji,jj) * r1_hu(ji,jj,Kmm_a)  & 
     996                      & + ( hdivdiff(ji+1,jj) - ze1v ) * r1_e1u(ji,jj)                       &  
     997                      & - rn_trelax_dyn * r1_Dt * fspu_2d(ji,jj) * ubdiff(ji,jj) 
     998 
     999                  ! add it to the general momentum trends 
     1000                  uu(ji,jj,:,Krhs_a) = uu(ji,jj,:,Krhs_a) + zua                                  
     1001               ENDIF 
     1002            END DO 
     1003         END DO 
     1004 
     1005         tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE. 
     1006 
     1007         jmax = j2-1 
     1008         ind1 = jpjglo - ( nn_hls + nbghostcells + 2 )   ! North 
     1009         DO jj = mj0(ind1), mj1(ind1)                  
     1010            jmax = MIN(jmax,jj) 
     1011         END DO 
     1012 
     1013         DO jj = j1+1, jmax 
     1014            DO ji = i1+1, i2   ! vector opt. 
     1015               IF (.NOT. tabspongedone_v(ji,jj)) THEN 
     1016                     ze2u = rotdiff (ji,jj) 
     1017                     ze1v = hdivdiff(ji,jj) 
     1018                     zva = + ( ze2u - rotdiff (ji-1,jj) ) * r1_e1v(ji,jj) * r1_hv(ji,jj,Kmm_a) & 
     1019                           + ( hdivdiff(ji,jj+1) - ze1v ) * r1_e2v(ji,jj) 
     1020                     vv(ji,jj,:,Krhs_a) = vv(ji,jj,:,Krhs_a) + zva 
     1021               ENDIF 
     1022            END DO 
     1023         END DO 
     1024         ! 
     1025         tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE. 
     1026         ! 
     1027      ENDIF 
     1028      ! 
     1029   END SUBROUTINE interpunb_sponge 
     1030 
     1031    
     1032   SUBROUTINE interpvnb_sponge(tabres,i1,i2,j1,j2, before) 
     1033      !!--------------------------------------------- 
     1034      !!   *** ROUTINE interpvnb_sponge *** 
     1035      !!---------------------------------------------  
     1036      INTEGER, INTENT(in) :: i1,i2,j1,j2 
     1037      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres 
     1038      LOGICAL, INTENT(in) :: before 
     1039      ! 
     1040      INTEGER  ::   ji, jj, ind1, imax 
     1041      REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr 
     1042      REAL(wp), DIMENSION(i1:i2,j1:j2) :: vbdiff 
     1043      REAL(wp), DIMENSION(i1:i2,j1:j2) :: rotdiff, hdivdiff 
     1044      !!---------------------------------------------  
     1045       
     1046      IF( before ) THEN  
     1047         DO jj=j1,j2 
     1048            DO ji=i1,i2 
     1049               tabres(ji,jj) = vv_b(ji,jj,Kmm_a) 
     1050            END DO 
     1051         END DO 
     1052      ELSE 
     1053         vbdiff(i1:i2,j1:j2) = (vv_b(i1:i2,j1:j2,Kmm_a) - tabres(i1:i2,j1:j2))*vmask(i1:i2,j1:j2,1) 
     1054         !                                             ! -------- 
     1055         ! Horizontal divergence                       !   div 
     1056         !                                             ! -------- 
     1057         DO jj = j1+1,j2 
     1058            DO ji = i1,i2   ! vector opt. 
     1059               zbtr = rn_sponge_dyn * r1_Dt * fspt_2d(ji,jj) * r1_ht_0(ji,jj) 
     1060               hdivdiff(ji,jj) = ( e1v(ji,jj  ) * hv(ji,jj  ,Kbb_a) * vbdiff(ji,jj  )  & 
     1061                               &  -e1v(ji,jj-1) * hv(ji,jj-1,Kbb_a) * vbdiff(ji,jj-1)  ) * zbtr 
     1062            END DO 
     1063         END DO 
     1064         DO jj = j1,j2 
     1065            DO ji = i1,i2-1   ! vector opt. 
     1066               zbtr = rn_sponge_dyn * r1_Dt * fspf_2d(ji,jj) * hf_0(ji,jj)  
     1067               rotdiff(ji,jj) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj) &  
     1068                              &  -e2v(ji  ,jj) * vbdiff(ji  ,jj)  ) * fmask(ji,jj,1) * zbtr 
     1069            END DO 
     1070         END DO 
     1071         !                                                ! =============== 
     1072         !                                                 
     1073 
     1074         imax = i2 - 1 
     1075         ind1 = jpiglo - ( nn_hls + nbghostcells + 2 )   ! East 
     1076         DO ji = mi0(ind1), mi1(ind1)                 
     1077            imax = MIN(imax,ji) 
     1078         END DO 
     1079          
     1080         DO jj = j1+1, j2 
     1081            DO ji = i1+1, imax   ! vector opt. 
     1082               IF( .NOT. tabspongedone_u(ji,jj) ) THEN                                                      
     1083                  zua = - ( rotdiff (ji  ,jj) - rotdiff (ji,jj-1)) * r1_e2u(ji,jj) * r1_hu(ji,jj,Kmm_a)  & 
     1084                      & + ( hdivdiff(ji+1,jj) - hdivdiff(ji,jj  )) * r1_e1u(ji,jj) 
     1085                  uu(ji,jj,:,Krhs_a) = uu(ji,jj,:,Krhs_a) + zua 
     1086               ENDIF 
     1087            END DO 
     1088         END DO 
     1089         ! 
     1090         tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE. 
     1091         ! 
     1092         DO jj = j1+1, j2-1 
     1093            DO ji = i1+1, i2-1   ! vector opt. 
     1094               IF( .NOT. tabspongedone_v(ji,jj) ) THEN 
     1095                  zva  =  ( rotdiff (ji,jj  ) - rotdiff (ji-1,jj) ) * r1_e1v(ji,jj) *r1_hv(ji,jj,Kmm_a) & 
     1096                     &  + ( hdivdiff(ji,jj+1) - hdivdiff(ji  ,jj) ) * r1_e2v(ji,jj)                     & 
     1097                     &  - rn_trelax_dyn * r1_Dt * fspv_2d(ji,jj) * vbdiff(ji,jj) 
     1098                  vv(ji,jj,:,Krhs_a) = vv(ji,jj,:,Krhs_a) + zva 
     1099               ENDIF 
     1100            END DO 
     1101         END DO 
     1102         tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE. 
     1103      ENDIF 
     1104      ! 
     1105   END SUBROUTINE interpvnb_sponge 
     1106 
    7881107 
    7891108#else 
  • NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_oce_update.F90

    r13674 r13937  
    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*/ 
     
    5151      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update tracers  from grid Number',Agrif_Fixed() 
    5252 
    53       Agrif_UseSpecialValueInUpdate = .NOT.ln_vert_remap 
     53      l_vremap                      = ln_vert_remap 
     54      Agrif_UseSpecialValueInUpdate = .NOT.l_vremap 
    5455      Agrif_SpecialValueFineGrid    = 0._wp 
    55       l_vremap                      = ln_vert_remap 
    5656      !  
    5757# if ! defined DECAL_FEEDBACK 
     
    8484      l_vremap                      = ln_vert_remap 
    8585      use_sign_north                = .TRUE. 
    86       sign_north                    = -1._wp 
    87  
    88       !      
     86      sign_north                    = -1._wp      
     87! 
     88# if ! defined DECAL_FEEDBACK_2D 
     89      CALL Agrif_Update_Variable(unb_update_id,locupdate1=(/  nn_shift_bar,-2/),locupdate2=(/  nn_shift_bar,-2/),procname = updateU2d) 
     90      CALL Agrif_Update_Variable(vnb_update_id,locupdate1=(/  nn_shift_bar,-2/),locupdate2=(/  nn_shift_bar,-2/),procname = updateV2d) 
     91# else 
     92      CALL Agrif_Update_Variable(unb_update_id,locupdate1=(/  nn_shift_bar,-2/),locupdate2=(/1+nn_shift_bar,-2/),procname = updateU2d) 
     93      CALL Agrif_Update_Variable(vnb_update_id,locupdate1=(/1+nn_shift_bar,-2/),locupdate2=(/  nn_shift_bar,-2/),procname = updateV2d)   
     94# endif 
     95      !  
     96      IF ( ln_dynspg_ts .AND. ln_bt_fw ) THEN 
     97         ! Update time integrated transports 
     98#  if ! defined DECAL_FEEDBACK_2D 
     99         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/  nn_shift_bar,-2/),locupdate2=(/  nn_shift_bar,-2/),procname = updateub2b) 
     100         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/  nn_shift_bar,-2/),locupdate2=(/  nn_shift_bar,-2/),procname = updatevb2b) 
     101#  else 
     102         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/  nn_shift_bar,-2/),locupdate2=(/1+nn_shift_bar,-2/),procname = updateub2b) 
     103         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1+nn_shift_bar,-2/),locupdate2=(/  nn_shift_bar,-2/),procname = updatevb2b) 
     104#  endif 
     105      END IF 
     106 
    89107# if ! defined DECAL_FEEDBACK 
    90108      CALL Agrif_Update_Variable(un_update_id,procname = updateU) 
     
    100118!      CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updateV) 
    101119# endif 
    102  
    103 # if ! defined DECAL_FEEDBACK_2D 
    104       CALL Agrif_Update_Variable(e1u_id,procname = updateU2d) 
    105       CALL Agrif_Update_Variable(e2v_id,procname = updateV2d)   
    106 # else 
    107       CALL Agrif_Update_Variable(e1u_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU2d) 
    108       CALL Agrif_Update_Variable(e2v_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV2d)   
    109 # endif 
    110       ! 
    111 # if ! defined DECAL_FEEDBACK_2D 
    112       ! Account for updated thicknesses at boundary edges 
    113       IF (.NOT.ln_linssh) THEN 
    114 !         CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_u_bdy) 
    115 !         CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_v_bdy) 
    116       ENDIF 
    117 # endif 
    118       !  
    119       IF ( ln_dynspg_ts .AND. ln_bt_fw ) THEN 
    120          ! Update time integrated transports 
    121 #  if ! defined DECAL_FEEDBACK_2D 
    122          CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b) 
    123          CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b) 
    124 #  else 
    125          CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateub2b) 
    126          CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updatevb2b) 
    127 #  endif 
    128       END IF 
    129120      ! 
    130121      use_sign_north = .FALSE. 
     
    143134      Agrif_SpecialValueFineGrid = 0._wp 
    144135# if ! defined DECAL_FEEDBACK_2D 
    145       CALL Agrif_Update_Variable(sshn_id,procname = updateSSH) 
     136      CALL Agrif_Update_Variable(sshn_id,locupdate=(/  nn_shift_bar,-2/), procname = updateSSH)  
    146137# else 
    147       CALL Agrif_Update_Variable(sshn_id,locupdate=(/1,0/),procname = updateSSH) 
     138      CALL Agrif_Update_Variable(sshn_id,locupdate=(/1+nn_shift_bar,-2/),procname = updateSSH) 
    148139# endif 
    149140      ! 
     
    151142      ! 
    152143#  if defined VOL_REFLUX 
    153       IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN 
     144      IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN  
    154145         use_sign_north = .TRUE. 
    155146         sign_north = -1._wp 
    156147         ! Refluxing on ssh: 
    157148#  if defined DECAL_FEEDBACK_2D 
    158          CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0, 0/),locupdate2=(/1, 1/),procname = reflux_sshu) 
    159          CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1, 1/),locupdate2=(/0, 0/),procname = reflux_sshv) 
     149         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) 
     150         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) 
    160151#  else 
    161          CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/-1,-1/),locupdate2=(/ 0, 0/),procname = reflux_sshu) 
    162          CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/ 0, 0/),locupdate2=(/-1,-1/),procname = reflux_sshv) 
     152         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) 
     153         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) 
    163154#  endif 
    164155         use_sign_north = .FALSE. 
     
    548539         END DO 
    549540         ! 
     541         ! Correct now and before transports: 
     542         DO jj=j1,j2 
     543            DO ji=i1,i2 
     544               spgu(ji,jj) = 0._wp 
     545               DO jk=1,jpkm1 
     546                  spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a) 
     547               END DO 
     548               ! 
     549               DO jk=1,jpkm1               
     550                  uu(ji,jj,jk,Kmm_a) = uu(ji,jj,jk,Kmm_a) + & 
     551                      &  (uu_b(ji,jj,Kmm_a) - spgu(ji,jj) * r1_hu(ji,jj,Kmm_a)) * umask(ji,jj,jk)            
     552               END DO 
     553               ! 
     554               spgu(ji,jj) = 0._wp 
     555               DO jk=1,jpkm1 
     556                  spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kbb_a) * uu(ji,jj,jk,Kbb_a) 
     557               END DO 
     558               ! 
     559               DO jk=1,jpkm1               
     560                  uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) + & 
     561                      &  (uu_b(ji,jj,Kbb_a) - spgu(ji,jj) * r1_hu(ji,jj,Kbb_a)) * umask(ji,jj,jk)            
     562               END DO 
     563               ! 
     564            END DO 
     565         END DO 
     566         ! 
    550567         IF  ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 
    551568            uu(i1:i2,j1:j2,1:jpkm1,Kbb_a)  = uu(i1:i2,j1:j2,1:jpkm1,Kmm_a) 
     
    555572      !  
    556573   END SUBROUTINE updateu 
    557  
    558    SUBROUTINE correct_u_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) 
    559       !!--------------------------------------------- 
    560       !!           *** ROUTINE correct_u_bdy *** 
    561       !!--------------------------------------------- 
    562       INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
    563       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    564       LOGICAL                                     , INTENT(in   ) :: before 
    565       INTEGER                                     , INTENT(in)    :: nb, ndir 
    566       !! 
    567       LOGICAL :: western_side, eastern_side  
    568       ! 
    569       INTEGER  :: jj, jk 
    570       REAL(wp) :: zcor 
    571       !!--------------------------------------------- 
    572       !  
    573       IF( .NOT.before ) THEN 
    574          ! 
    575          western_side  = (nb == 1).AND.(ndir == 1) 
    576          eastern_side  = (nb == 1).AND.(ndir == 2) 
    577          ! 
    578          IF (western_side) THEN 
    579             DO jj=j1,j2 
    580                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) 
    581                uu_b(i1-1,jj,Kmm_a) = uu_b(i1-1,jj,Kmm_a) + zcor 
    582                DO jk=1,jpkm1 
    583                   uu(i1-1,jj,jk,Kmm_a) = uu(i1-1,jj,jk,Kmm_a) + zcor * umask(i1-1,jj,jk) 
    584                END DO  
    585             END DO 
    586          ENDIF 
    587          ! 
    588          IF (eastern_side) THEN 
    589             DO jj=j1,j2 
    590                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) 
    591                uu_b(i2+1,jj,Kmm_a) = uu_b(i2+1,jj,Kmm_a) + zcor 
    592                DO jk=1,jpkm1 
    593                   uu(i2+1,jj,jk,Kmm_a) = uu(i2+1,jj,jk,Kmm_a) + zcor * umask(i2+1,jj,jk) 
    594                END DO  
    595             END DO 
    596          ENDIF 
    597          ! 
    598       ENDIF 
    599       !  
    600    END SUBROUTINE correct_u_bdy 
    601574 
    602575 
     
    712685         END DO 
    713686         ! 
     687         ! Correct now and before transports: 
     688         DO jj=j1,j2 
     689            DO ji=i1,i2 
     690               spgv(ji,jj) = 0._wp 
     691               DO jk=1,jpkm1 
     692                  spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a) 
     693               END DO 
     694               ! 
     695               DO jk=1,jpkm1               
     696                  vv(ji,jj,jk,Kmm_a) = vv(ji,jj,jk,Kmm_a) + & 
     697                      &  (vv_b(ji,jj,Kmm_a) - spgv(ji,jj) * r1_hv(ji,jj,Kmm_a)) * vmask(ji,jj,jk)            
     698               END DO 
     699               ! 
     700               spgv(ji,jj) = 0._wp 
     701               DO jk=1,jpkm1 
     702                  spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kbb_a) * vv(ji,jj,jk,Kbb_a) 
     703               END DO 
     704               ! 
     705               DO jk=1,jpkm1               
     706                  vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) + & 
     707                      &  (vv_b(ji,jj,Kbb_a) - spgv(ji,jj) * r1_hv(ji,jj,Kbb_a)) * vmask(ji,jj,jk)            
     708               END DO 
     709               ! 
     710            END DO 
     711         END DO 
     712         ! 
    714713         IF  ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 
    715714            vv(i1:i2,j1:j2,1:jpkm1,Kbb_a)  = vv(i1:i2,j1:j2,1:jpkm1,Kmm_a) 
     
    719718      !  
    720719   END SUBROUTINE updatev 
    721  
    722  
    723    SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) 
    724       !!--------------------------------------------- 
    725       !!           *** ROUTINE correct_v_bdy *** 
    726       !!--------------------------------------------- 
    727       INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
    728       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    729       LOGICAL                                     , INTENT(in   ) :: before 
    730       INTEGER                                     , INTENT(in)    :: nb, ndir 
    731       !! 
    732       LOGICAL :: southern_side, northern_side  
    733       ! 
    734       INTEGER  :: ji, jk 
    735       REAL(wp) :: zcor 
    736       !!--------------------------------------------- 
    737       !  
    738       IF( .NOT.before ) THEN 
    739          ! 
    740          southern_side = (nb == 2).AND.(ndir == 1) 
    741          northern_side = (nb == 2).AND.(ndir == 2) 
    742          ! 
    743          IF (southern_side) THEN 
    744             DO ji=i1,i2 
    745                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) 
    746                vv_b(ji,j1-1,Kmm_a) = vv_b(ji,j1-1,Kmm_a) + zcor 
    747                DO jk=1,jpkm1 
    748                   vv(ji,j1-1,jk,Kmm_a) = vv(ji,j1-1,jk,Kmm_a) + zcor * vmask(ji,j1-1,jk) 
    749                END DO  
    750             END DO 
    751          ENDIF 
    752          ! 
    753          IF (northern_side) THEN 
    754             DO ji=i1,i2 
    755                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) 
    756                vv_b(ji,j2+1,Kmm_a) = vv_b(ji,j2+1,Kmm_a) + zcor 
    757                DO jk=1,jpkm1 
    758                   vv(ji,j2+1,jk,Kmm_a) = vv(ji,j2+1,jk,Kmm_a) + zcor * vmask(ji,j2+1,jk) 
    759                END DO  
    760             END DO 
    761          ENDIF 
    762          ! 
    763       ENDIF 
    764       !  
    765    END SUBROUTINE correct_v_bdy 
    766720 
    767721 
     
    791745               tabres(ji,jj) =  tabres(ji,jj) * r1_e2u(ji,jj)   
    792746               !     
    793                ! Update "now" 3d velocities: 
    794                spgu(ji,jj) = 0._wp 
    795                DO jk=1,jpkm1 
    796                   spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a) 
    797                END DO 
    798                ! 
    799                zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu(ji,jj,Kmm_a) 
    800                DO jk=1,jpkm1               
    801                   uu(ji,jj,jk,Kmm_a) = uu(ji,jj,jk,Kmm_a) + zcorr * umask(ji,jj,jk)            
    802                END DO 
    803                ! 
    804747               ! Update barotropic velocities: 
    805748               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 
     
    809752                  END IF 
    810753               ENDIF     
    811                uu_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hu(ji,jj,Kmm_a) * umask(ji,jj,1) 
    812                !        
    813                ! Correct "before" velocities to hold correct bt component: 
    814                spgu(ji,jj) = 0.e0 
    815                DO jk=1,jpkm1 
    816                   spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kbb_a) * uu(ji,jj,jk,Kbb_a) 
    817                END DO 
    818                ! 
    819                zcorr = uu_b(ji,jj,Kbb_a) - spgu(ji,jj) * r1_hu(ji,jj,Kbb_a) 
    820                DO jk=1,jpkm1               
    821                   uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) + zcorr * umask(ji,jj,jk)            
    822                END DO 
     754               uu_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hu(ji,jj,Kmm_a) * umask(ji,jj,1)        
    823755               ! 
    824756            END DO 
     
    855787         DO jj=j1,j2 
    856788            DO ji=i1,i2 
    857                tabres(ji,jj) =  tabres(ji,jj) * r1_e1v(ji,jj)   
    858                !     
    859                ! Update "now" 3d velocities: 
    860                spgv(ji,jj) = 0.e0 
    861                DO jk=1,jpkm1 
    862                   spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a) 
    863                END DO 
    864                ! 
    865                zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv(ji,jj,Kmm_a) 
    866                DO jk=1,jpkm1               
    867                   vv(ji,jj,jk,Kmm_a) = vv(ji,jj,jk,Kmm_a) + zcorr * vmask(ji,jj,jk)            
    868                END DO 
     789               tabres(ji,jj) =  tabres(ji,jj) * r1_e1v(ji,jj)       
    869790               ! 
    870791               ! Update barotropic velocities: 
     
    875796                  END IF 
    876797               ENDIF               
    877                vv_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hv(ji,jj,Kmm_a) * vmask(ji,jj,1) 
    878                !        
    879                ! Correct "before" velocities to hold correct bt component: 
    880                spgv(ji,jj) = 0.e0 
    881                DO jk=1,jpkm1 
    882                   spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kbb_a) * vv(ji,jj,jk,Kbb_a) 
    883                END DO 
    884                ! 
    885                zcorr = vv_b(ji,jj,Kbb_a) - spgv(ji,jj) * r1_hv(ji,jj,Kbb_a) 
    886                DO jk=1,jpkm1               
    887                   vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) + zcorr * vmask(ji,jj,jk)            
    888                END DO 
     798               vv_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hv(ji,jj,Kmm_a) * vmask(ji,jj,1)      
    889799               ! 
    890800            END DO 
     
    13431253      IF(lwp) WRITE(numout,*) 'AGRIF: Check parent volume at Level:', Agrif_Level() 
    13441254      ! 
    1345 # if ! defined DECAL_FEEDBACK && ! defined DECAL_FEEDBACK_2D 
     1255# if ! defined DECAL_FEEDBACK 
    13461256      CALL Agrif_Update_Variable(batupd_id,procname = update_bat) 
    13471257# else 
  • NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_user.F90

    r13673 r13937  
    11#undef UPD_HIGH   /* MIX HIGH UPDATE */ 
     2#define DIV_CONS   /* DIVERGENCE CONS */ 
    23#if defined key_agrif 
    34   !! * Substitutions 
     
    5556      IMPLICIT NONE 
    5657      ! 
    57       INTEGER :: ind1, ind2, ind3 
     58      INTEGER :: ind1, ind2, ind3, imaxrho 
    5859      INTEGER :: its 
    5960      External :: nemo_mapping 
     
    7778      ! 1. Declaration of the type of variable which have to be interpolated 
    7879      !--------------------------------------------------------------------- 
    79       ind1 =              nbghostcells 
     80      ind1 =              nbghostcells  
    8081      ind2 = nn_hls + 2 + nbghostcells_x 
    8182      ind3 = nn_hls + 2 + nbghostcells_y_s 
    82  
    83       CALL agrif_declare_variable((/2,2,0  /),(/ind2  ,ind3,0    /),(/'x','y','N'    /),(/1,1,1  /),(/jpi,jpj,jpk    /),   e3t_id) 
    84       CALL agrif_declare_variable((/2,2    /),(/ind2  ,ind3      /),(/'x','y'        /),(/1,1    /),(/jpi,jpj        /),  mbkt_id) 
    85       CALL agrif_declare_variable((/2,2    /),(/ind2  ,ind3      /),(/'x','y'        /),(/1,1    /),(/jpi,jpj        /),   ht0_id) 
    86  
    87       CALL agrif_declare_variable((/1,2    /),(/ind2-1,ind3      /),(/'x','y'        /),(/1,1    /),(/jpi,jpj        /),   e1u_id) 
    88       CALL agrif_declare_variable((/2,1    /),(/ind2  ,ind3-1    /),(/'x','y'        /),(/1,1    /),(/jpi,jpj        /),   e2v_id) 
     83      imaxrho = MAX(Agrif_irhox(), Agrif_irhoy()) 
     84 
     85      CALL agrif_declare_variable((/2,2,0  /),(/ind2  ,ind3,0    /),(/'x','y','N'    /),(/1,1,1  /),(/jpi,jpj,jpk    /),        e3t_id) 
     86      CALL agrif_declare_variable((/2,2,0  /),(/ind2  ,ind3,0    /),(/'x','y','N'    /),(/1,1,1  /),(/jpi,jpj,jpk    /),e3t0_interp_id) 
     87      CALL agrif_declare_variable((/2,2    /),(/ind2  ,ind3      /),(/'x','y'        /),(/1,1    /),(/jpi,jpj        /),       mbkt_id) 
     88      CALL agrif_declare_variable((/2,2    /),(/ind2  ,ind3      /),(/'x','y'        /),(/1,1    /),(/jpi,jpj        /),        ht0_id) 
    8989    
    9090      ! Initial or restart velues 
     
    100100      ! 2. Type of interpolation 
    101101      !------------------------- 
    102       CALL Agrif_Set_bcinterp(   e3t_id,interp =AGRIF_constant) 
    103  
    104       CALL Agrif_Set_bcinterp(  mbkt_id,interp =AGRIF_constant) 
    105       CALL Agrif_Set_interp  (  mbkt_id,interp =AGRIF_constant) 
    106       CALL Agrif_Set_bcinterp(   ht0_id,interp =AGRIF_constant) 
    107       CALL Agrif_Set_interp  (   ht0_id,interp =AGRIF_constant) 
    108  
    109       CALL Agrif_Set_bcinterp(   e1u_id,interp1=Agrif_linear, interp2=AGRIF_ppm    ) 
    110       CALL Agrif_Set_bcinterp(   e2v_id,interp1=AGRIF_ppm   , interp2=Agrif_linear ) 
     102      CALL Agrif_Set_bcinterp(        e3t_id,interp =AGRIF_constant) 
     103      CALL Agrif_Set_bcinterp(e3t0_interp_id,interp =AGRIF_linear  ) 
     104      CALL Agrif_Set_interp  (e3t0_interp_id,interp =AGRIF_linear  ) 
     105      CALL Agrif_Set_bcinterp(       mbkt_id,interp =AGRIF_constant) 
     106      CALL Agrif_Set_interp  (       mbkt_id,interp =AGRIF_constant) 
     107      CALL Agrif_Set_bcinterp(        ht0_id,interp =AGRIF_constant) 
     108      CALL Agrif_Set_interp  (        ht0_id,interp =AGRIF_constant) 
    111109 
    112110      ! Initial fields 
     
    122120       ! 3. Location of interpolation 
    123121      !----------------------------- 
    124 !      CALL Agrif_Set_bc(  e3t_id, (/-nn_sponge_len*Agrif_irhox(),ind1-1/) )   
     122!      CALL Agrif_Set_bc(  e3t_id, (/-nn_sponge_len*imaxrho,ind1-1/) )   
    125123! JC: check near the boundary only until matching in sponge has been sorted out: 
    126124      CALL Agrif_Set_bc(    e3t_id, (/0,ind1-1/) )   
     
    128126      ! extend the interpolation zone by 1 more point than necessary: 
    129127      ! RB check here 
    130       CALL Agrif_Set_bc(   mbkt_id, (/-nn_sponge_len*Agrif_irhox()-2,ind1/) ) 
    131       CALL Agrif_Set_bc(    ht0_id, (/-nn_sponge_len*Agrif_irhox()-2,ind1/) ) 
    132        
    133       CALL Agrif_Set_bc(    e1u_id, (/0,ind1-1/) ) 
    134       CALL Agrif_Set_bc(    e2v_id, (/0,ind1-1/) )   
    135  
    136       CALL Agrif_Set_bc(  tsini_id, (/0,ind1-1/) ) ! if west,  rhox=3 and nbghost=3: columns 2 to 4 
    137       CALL Agrif_Set_bc(   uini_id, (/0,ind1-1/) )  
    138       CALL Agrif_Set_bc(   vini_id, (/0,ind1-1/) ) 
    139       CALL Agrif_Set_bc( sshini_id, (/0,ind1-1/) ) 
     128      CALL Agrif_Set_bc( e3t0_interp_id, (/-nn_sponge_len*imaxrho-2,ind1/) ) 
     129      CALL Agrif_Set_bc(        mbkt_id, (/-nn_sponge_len*imaxrho-2,ind1/) ) 
     130      CALL Agrif_Set_bc(         ht0_id, (/-nn_sponge_len*imaxrho-2,ind1/) ) 
     131      CALL Agrif_Set_bc(       tsini_id, (/0,ind1-1/) ) ! if west,  rhox=3 and nbghost=3: columns 2 to 4 
     132      CALL Agrif_Set_bc(        uini_id, (/0,ind1-1/) )  
     133      CALL Agrif_Set_bc(        vini_id, (/0,ind1-1/) ) 
     134      CALL Agrif_Set_bc(      sshini_id, (/0,ind1-1/) ) 
    140135 
    141136      ! 4. Update type 
    142137      !---------------  
    143138# if defined UPD_HIGH 
    144       CALL Agrif_Set_Updatetype(e1u_id,update1 = Agrif_Update_Average       , update2=Agrif_Update_Full_Weighting) 
    145       CALL Agrif_Set_Updatetype(e2v_id,update1 = Agrif_Update_Full_Weighting, update2=Agrif_Update_Average       ) 
    146139      CALL Agrif_Set_Updatetype(batupd_id, update = Agrif_Update_Full_Weighting) 
    147140#else 
    148       CALL Agrif_Set_Updatetype(e1u_id,update1 = Agrif_Update_Copy          , update2=Agrif_Update_Average       ) 
    149       CALL Agrif_Set_Updatetype(e2v_id,update1 = Agrif_Update_Average       , update2=Agrif_Update_Copy          ) 
    150141      CALL Agrif_Set_Updatetype(batupd_id, update = Agrif_Update_Average) 
    151142#endif       
     
    181172      REAL(wp), DIMENSION(jpi,jpj) ::   zk   ! workspace 
    182173      INTEGER :: ji, jj, jk 
     174      INTEGER :: jpk_parent, ierr 
    183175      !!---------------------------------------------------------------------- 
    184176     
     
    227219      END_2D 
    228220      CALL lbc_lnk( 'Agrif_InitValues_Domain', zk, 'V', 1.0_wp ) 
    229       mbkv_parent(:,:) = MAX( NINT( zk(:,:) ), 1 )    
    230  
     221      mbkv_parent(:,:) = MAX( NINT( zk(:,:) ), 1 )   
     222      ! 
     223      ! Build "intermediate" parent vertical grid on child domain 
     224      IF ( ln_vert_remap ) THEN 
     225 
     226         jpk_parent = Agrif_parent( jpk ) 
     227         ALLOCATE(e3t0_parent(jpi,jpj,jpk_parent), & 
     228            &     e3u0_parent(jpi,jpj,jpk_parent), & 
     229            &     e3v0_parent(jpi,jpj,jpk_parent), STAT = ierr)  
     230         IF( ierr  > 0 )   CALL ctl_warn('Agrif_Init_Domain: allocation of arrays failed') 
     231        
     232         ! Retrieve expected parent scale factors on child grid: 
     233         Agrif_UseSpecialValue = .FALSE. 
     234         e3t0_parent(:,:,:) = 0._wp 
     235         CALL Agrif_Init_Variable(e3t0_interp_id, procname=interpe3t0_vremap) 
     236         ! 
     237         ! Deduce scale factors at U and V points: 
     238         DO_3D( 0, 0, 0, 0, 1, jpk_parent ) 
     239            e3u0_parent(ji,jj,jk) = 0.5_wp * (e3t0_parent(ji,jj,jk) + e3t0_parent(ji+1,jj  ,jk)) 
     240            e3v0_parent(ji,jj,jk) = 0.5_wp * (e3t0_parent(ji,jj,jk) + e3t0_parent(ji  ,jj+1,jk)) 
     241         END_3D 
     242 
     243         ! Assume a step at the bottom except if (pure) s-coordinates 
     244         IF ( .NOT.Agrif_Parent(ln_sco) ) THEN  
     245            DO_2D( 1, 0, 1, 0 ) 
     246               jk = mbku_parent(ji,jj) 
     247               e3u0_parent(ji,jj,jk) = MIN(e3t0_parent(ji,jj,jk), e3t0_parent(ji+1,jj  ,jk)) 
     248               jk = mbkv_parent(ji,jj) 
     249               e3v0_parent(ji,jj,jk) = MIN(e3t0_parent(ji,jj,jk), e3t0_parent(ji  ,jj+1,jk)) 
     250            END_2D 
     251         ENDIF 
     252 
     253         CALL lbc_lnk_multi( 'Agrif_Init_Domain', e3u0_parent, 'U', 1.0_wp, e3v0_parent, 'V', 1.0_wp ) 
     254      ENDIF 
    231255 
    232256      ! check if masks and bathymetries match 
     
    259283      ENDIF 
    260284      ! 
     285      WHERE (ssumask(:,:) == 0._wp) mbku_parent(:,:) = 0 
     286      WHERE (ssvmask(:,:) == 0._wp) mbkv_parent(:,:) = 0 
     287      WHERE (ssmask(:,:)  == 0._wp) mbkt_parent(:,:) = 0 
     288      ! 
    261289   END SUBROUTINE Agrif_Init_Domain 
    262290 
     
    315343      tabspongedone_v = .FALSE. 
    316344      CALL Agrif_Bc_variable(vn_sponge_id,calledweight=1.,procname=interpvn_sponge) 
     345      IF (nn_shift_bar>0) THEN 
     346         CALL Agrif_Sponge_2d 
     347         tabspongedone_u = .FALSE. 
     348         tabspongedone_v = .FALSE. 
     349         CALL Agrif_Bc_variable(unb_sponge_id,calledweight=1.,procname=interpunb_sponge) 
     350         tabspongedone_u = .FALSE. 
     351         tabspongedone_v = .FALSE. 
     352         CALL Agrif_Bc_variable(vnb_sponge_id,calledweight=1.,procname=interpvnb_sponge) 
     353      ENDIF 
    317354      use_sign_north = .FALSE. 
    318355      uu(:,:,:,Krhs_a) = 0._wp 
     
    328365         use_sign_north = .TRUE. 
    329366         sign_north = -1. 
    330          CALL Agrif_Bc_variable(        unb_id,calledweight=1.,procname=interpunb ) 
    331          CALL Agrif_Bc_variable(        vnb_id,calledweight=1.,procname=interpvnb ) 
     367         CALL Agrif_Bc_variable( unb_interp_id,calledweight=1.,procname=interpunb ) 
     368         CALL Agrif_Bc_variable( vnb_interp_id,calledweight=1.,procname=interpvnb ) 
    332369         CALL Agrif_Bc_variable(ub2b_interp_id,calledweight=1.,procname=interpub2b) 
    333370         CALL Agrif_Bc_variable(vb2b_interp_id,calledweight=1.,procname=interpvb2b) 
     371         use_sign_north = .FALSE. 
     372         ubdy(:,:) = 0._wp 
     373         vbdy(:,:) = 0._wp 
     374      ELSE 
     375         Agrif_UseSpecialValue = ln_spc_dyn 
     376         use_sign_north = .TRUE. 
     377         sign_north = -1. 
     378         CALL Agrif_Bc_variable( unb_interp_id,calledweight=1.,procname=interpunb ) 
     379         CALL Agrif_Bc_variable( vnb_interp_id,calledweight=1.,procname=interpvnb ) 
    334380         use_sign_north = .FALSE. 
    335381         ubdy(:,:) = 0._wp 
     
    386432      IMPLICIT NONE 
    387433      ! 
    388       INTEGER :: ind1, ind2, ind3 
     434      INTEGER :: ind1, ind2, ind3, imaxrho 
    389435      !!---------------------------------------------------------------------- 
    390436 
     
    394440      ind2 = nn_hls + 2 + nbghostcells_x 
    395441      ind3 = nn_hls + 2 + nbghostcells_y_s 
     442      imaxrho = MAX(Agrif_irhox(), Agrif_irhoy()) 
    396443 
    397444      CALL agrif_declare_variable((/2,2,0,0/),(/ind2,ind3,0,0/)  ,(/'x','y','N','N'/),(/1,1,1,1/),(/jpi,jpj,jpk,jpts+1/),ts_interp_id) 
     
    405452      CALL agrif_declare_variable((/2,1,0,0/),(/ind2,ind3-1,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/jpi,jpj,jpk,2/),vn_sponge_id) 
    406453 
    407       CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/jpi,jpj/),unb_id) 
    408       CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/jpi,jpj/),vnb_id) 
     454      CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/jpi,jpj/)  ,sshn_id) 
     455      CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/jpi,jpj/), unb_interp_id) 
     456      CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/jpi,jpj/), vnb_interp_id) 
    409457      CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/jpi,jpj/),ub2b_interp_id) 
    410458      CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/jpi,jpj/),vb2b_interp_id) 
     459      CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/jpi,jpj/), unb_sponge_id) 
     460      CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/jpi,jpj/), vnb_sponge_id) 
    411461      CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/jpi,jpj/),ub2b_update_id) 
    412462      CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/jpi,jpj/),vb2b_update_id) 
    413  
     463      CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/jpi,jpj/), unb_update_id) 
     464      CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/jpi,jpj/), vnb_update_id) 
     465      CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/jpi,jpj/)  ,ub2b_cor_id) 
     466      CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/jpi,jpj/)  ,vb2b_cor_id) 
    414467!      CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/jpi,jpj/),glamt_id) 
    415468!      CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/jpi,jpj/),gphit_id) 
    416       CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/jpi,jpj/),sshn_id) 
    417469 
    418470 
     
    426478      !------------------------- 
    427479      CALL Agrif_Set_bcinterp( ts_interp_id,interp =AGRIF_linear) 
    428       CALL Agrif_Set_bcinterp( un_interp_id,interp1=Agrif_linear,interp2=AGRIF_ppm   ) 
    429       CALL Agrif_Set_bcinterp( vn_interp_id,interp1=AGRIF_ppm   ,interp2=Agrif_linear) 
    430  
    431       CALL Agrif_Set_bcinterp(  ts_sponge_id,interp =AGRIF_linear) 
    432       CALL Agrif_Set_bcinterp(  un_sponge_id,interp1=Agrif_linear,interp2=AGRIF_ppm   ) 
    433       CALL Agrif_Set_bcinterp(  vn_sponge_id,interp1=AGRIF_ppm   ,interp2=Agrif_linear) 
    434  
    435       CALL Agrif_Set_bcinterp(       sshn_id,interp =AGRIF_linear) 
    436       CALL Agrif_Set_bcinterp(        unb_id,interp1=Agrif_linear,interp2=AGRIF_ppm   ) 
    437       CALL Agrif_Set_bcinterp(        vnb_id,interp1=AGRIF_ppm   ,interp2=Agrif_linear) 
    438       CALL Agrif_Set_bcinterp(ub2b_interp_id,interp1=Agrif_linear,interp2=AGRIF_ppm   ) 
    439       CALL Agrif_Set_bcinterp(vb2b_interp_id,interp1=AGRIF_ppm   ,interp2=Agrif_linear) 
    440 ! 
    441 ! > Divergence conserving alternative: 
    442 !      CALL Agrif_Set_bcinterp( ts_interp_id,interp =AGRIF_constant) 
    443 !      CALL Agrif_Set_bcinterp( un_interp_id,interp1=Agrif_linear,interp2=AGRIF_constant   ) 
    444 !      CALL Agrif_Set_bcinterp( vn_interp_id,interp1=AGRIF_constant   ,interp2=Agrif_linear) 
    445 ! 
    446 !      CALL Agrif_Set_bcinterp(  ts_sponge_id,interp =AGRIF_constant) 
    447 !      CALL Agrif_Set_bcinterp(  un_sponge_id,interp1=Agrif_linear,interp2=AGRIF_constant   ) 
    448 !      CALL Agrif_Set_bcinterp(  vn_sponge_id,interp1=AGRIF_constant   ,interp2=Agrif_linear) 
    449 ! 
    450 !      CALL Agrif_Set_bcinterp(sshn_id,interp=AGRIF_constant) 
    451 !      CALL Agrif_Set_bcinterp(unb_id,interp1=Agrif_linear,interp2=AGRIF_constant) 
    452 !      CALL Agrif_Set_bcinterp(vnb_id,interp1=AGRIF_constant,interp2=Agrif_linear) 
    453 !      CALL Agrif_Set_bcinterp(ub2b_interp_id,interp1=Agrif_linear,interp2=AGRIF_constant) 
    454 !      CALL Agrif_Set_bcinterp(vb2b_interp_id,interp1=AGRIF_constant,interp2=Agrif_linear) 
    455 !< 
     480      CALL Agrif_Set_bcinterp( ts_sponge_id,interp =AGRIF_linear) 
     481 
     482#if defined DIV_CONS 
     483      lk_tint2d_notinterp = .TRUE. 
     484      CALL Agrif_Set_bcinterp(sshn_id,interp=AGRIF_constant) 
     485      CALL Agrif_Set_bcinterp(ub2b_cor_id,interp=AGRIF_constant) 
     486      CALL Agrif_Set_bcinterp(vb2b_cor_id,interp=AGRIF_constant) 
     487      CALL Agrif_Set_bcinterp(ub2b_interp_id,interp1=Agrif_linear,interp2=AGRIF_linearconserv) 
     488      CALL Agrif_Set_bcinterp(vb2b_interp_id,interp1=AGRIF_linearconserv,interp2=Agrif_linear) 
     489#else 
     490      lk_tint2d_notinterp = .FALSE. 
     491      CALL Agrif_Set_bcinterp(sshn_id,interp =AGRIF_linear) 
     492      CALL Agrif_Set_bcinterp(ub2b_interp_id,interp1=Agrif_linear,interp2=AGRIF_ppm)  
     493      CALL Agrif_Set_bcinterp(vb2b_interp_id,interp1=AGRIF_ppm,interp2=Agrif_linear) 
     494#endif 
     495 
     496      CALL Agrif_Set_bcinterp(unb_interp_id,interp1=Agrif_linear,interp2=AGRIF_ppm   ) 
     497      CALL Agrif_Set_bcinterp(vnb_interp_id,interp1=AGRIF_ppm   ,interp2=Agrif_linear) 
     498      CALL Agrif_Set_bcinterp(unb_sponge_id,interp1=Agrif_linear,interp2=AGRIF_ppm   ) 
     499      CALL Agrif_Set_bcinterp(vnb_sponge_id,interp1=AGRIF_ppm   ,interp2=Agrif_linear) 
     500 
     501      CALL Agrif_Set_bcinterp(un_interp_id,interp1=Agrif_linear,interp2=AGRIF_ppm   ) 
     502      CALL Agrif_Set_bcinterp(vn_interp_id,interp1=AGRIF_ppm   ,interp2=Agrif_linear) 
     503 
     504      CALL Agrif_Set_bcinterp(un_sponge_id,interp1=Agrif_linear,interp2=AGRIF_ppm   ) 
     505      CALL Agrif_Set_bcinterp(vn_sponge_id,interp1=AGRIF_ppm   ,interp2=Agrif_linear) 
    456506 
    457507      IF( ln_zdftke.OR.ln_zdfgls )  CALL Agrif_Set_bcinterp( avm_id, interp=AGRIF_linear ) 
     
    463513      ! 3. Location of interpolation 
    464514      !----------------------------- 
    465       CALL Agrif_Set_bc( ts_interp_id, (/0,ind1-1/) ) ! if west,  rhox=3 and nbghost=3: columns 2 to 4 
    466       CALL Agrif_Set_bc( un_interp_id, (/0,ind1-1/) )  
    467       CALL Agrif_Set_bc( vn_interp_id, (/0,ind1-1/) ) 
    468  
    469       CALL Agrif_Set_bc(  ts_sponge_id, (/-nn_sponge_len*Agrif_irhox()-1,0/) )  ! if west,  rhox=3, nn_sponge_len=2  
    470       CALL Agrif_Set_bc(  un_sponge_id, (/-nn_sponge_len*Agrif_irhox()-1,0/) )  ! and nbghost=3:  
    471       CALL Agrif_Set_bc(  vn_sponge_id, (/-nn_sponge_len*Agrif_irhox()-1,0/) )  ! columns 4 to 11 
    472  
    473       CALL Agrif_Set_bc(        sshn_id, (/0,ind1-1/) ) 
    474       CALL Agrif_Set_bc(         unb_id, (/0,ind1-1/) ) 
    475       CALL Agrif_Set_bc(         vnb_id, (/0,ind1-1/) ) 
    476       CALL Agrif_Set_bc( ub2b_interp_id, (/0,ind1-1/) ) 
    477       CALL Agrif_Set_bc( vb2b_interp_id, (/0,ind1-1/) ) 
    478  
     515      CALL Agrif_Set_bc(  ts_interp_id, (/0,ind1-1/) ) ! if west,  rhox=3 and nbghost=3: columns 2 to 4 
     516      CALL Agrif_Set_bc(  un_interp_id, (/0,ind1-1/) )  
     517      CALL Agrif_Set_bc(  vn_interp_id, (/0,ind1-1/) ) 
     518 
     519      CALL Agrif_Set_bc(  ts_sponge_id, (/-nn_sponge_len*imaxrho-1,0/) )  ! if west,  rhox=3, nn_sponge_len=2  
     520      CALL Agrif_Set_bc(  un_sponge_id, (/-nn_sponge_len*imaxrho-1,0/) )  ! and nbghost=3:  
     521      CALL Agrif_Set_bc(  vn_sponge_id, (/-nn_sponge_len*imaxrho-1,0/) )  ! columns 4 to 11 
     522 
     523      CALL Agrif_Set_bc(       sshn_id, (/-imaxrho*nn_shift_bar,ind1-1/) ) 
     524      CALL Agrif_Set_bc( unb_interp_id, (/-imaxrho*nn_shift_bar,ind1-1/) ) 
     525      CALL Agrif_Set_bc( vnb_interp_id, (/-imaxrho*nn_shift_bar,ind1-1/) ) 
     526      CALL Agrif_Set_bc(ub2b_interp_id, (/-imaxrho*nn_shift_bar,ind1-1/) ) 
     527      CALL Agrif_Set_bc(vb2b_interp_id, (/-imaxrho*nn_shift_bar,ind1-1/) ) 
     528      CALL Agrif_Set_bc( unb_sponge_id, (/-(nn_sponge_len+nn_shift_bar)*imaxrho,-imaxrho*nn_shift_bar/) ) 
     529      CALL Agrif_Set_bc( vnb_sponge_id, (/-(nn_sponge_len+nn_shift_bar)*imaxrho,-imaxrho*nn_shift_bar/) ) 
     530      CALL Agrif_Set_bc(   ub2b_cor_id, (/-imaxrho*nn_shift_bar,ind1/) ) 
     531      CALL Agrif_Set_bc(   vb2b_cor_id, (/-imaxrho*nn_shift_bar,ind1/) ) 
    479532      IF( ln_zdftke.OR.ln_zdfgls ) CALL Agrif_Set_bc( avm_id, (/0,ind1/) ) 
    480533!!$      CALL Agrif_Set_bc(glamt_id, (/0,ind1-1/) )   
     
    485538 
    486539# if defined UPD_HIGH 
    487       CALL Agrif_Set_Updatetype(ts_interp_id,update  = Agrif_Update_Full_Weighting) 
    488       CALL Agrif_Set_Updatetype(un_update_id,update1 = Agrif_Update_Average       , update2 = Agrif_Update_Full_Weighting) 
    489       CALL Agrif_Set_Updatetype(vn_update_id,update1 = Agrif_Update_Full_Weighting, update2 = Agrif_Update_Average       ) 
    490  
     540      CALL Agrif_Set_Updatetype(  ts_interp_id,update  = Agrif_Update_Full_Weighting) 
     541      CALL Agrif_Set_Updatetype(  un_update_id,update1 = Agrif_Update_Average       , update2 = Agrif_Update_Full_Weighting) 
     542      CALL Agrif_Set_Updatetype(  vn_update_id,update1 = Agrif_Update_Full_Weighting, update2 = Agrif_Update_Average       ) 
     543 
     544      CALL Agrif_Set_Updatetype( unb_update_id,update1 = Agrif_Update_Average       , update2 = Agrif_Update_Full_Weighting)  
     545      CALL Agrif_Set_Updatetype( vnb_update_id,update1 = Agrif_Update_Full_Weighting, update2 = Agrif_Update_Average       ) 
    491546      CALL Agrif_Set_Updatetype(ub2b_update_id,update1 = Agrif_Update_Average       , update2 = Agrif_Update_Full_Weighting) 
    492547      CALL Agrif_Set_Updatetype(vb2b_update_id,update1 = Agrif_Update_Full_Weighting, update2 = Agrif_Update_Average       ) 
     
    501556 
    502557#else 
    503       CALL Agrif_Set_Updatetype(ts_update_id ,update  = AGRIF_Update_Average) 
    504       CALL Agrif_Set_Updatetype(un_update_id ,update1 = Agrif_Update_Copy   , update2 = Agrif_Update_Average) 
    505       CALL Agrif_Set_Updatetype(vn_update_id ,update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy   ) 
    506  
     558      CALL Agrif_Set_Updatetype(  ts_update_id,update  = AGRIF_Update_Average) 
     559      CALL Agrif_Set_Updatetype(  un_update_id,update1 = Agrif_Update_Copy   , update2 = Agrif_Update_Average) 
     560      CALL Agrif_Set_Updatetype(  vn_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy   ) 
     561 
     562      CALL Agrif_Set_Updatetype( unb_update_id,update1 = Agrif_Update_Copy   , update2 = Agrif_Update_Average) 
     563      CALL Agrif_Set_Updatetype( vnb_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy   ) 
    507564      CALL Agrif_Set_Updatetype(ub2b_update_id,update1 = Agrif_Update_Copy   , update2 = Agrif_Update_Average) 
    508565      CALL Agrif_Set_Updatetype(vb2b_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy   ) 
     
    715772      IMPLICIT NONE 
    716773      ! 
    717       INTEGER :: ind1, ind2, ind3 
     774      INTEGER :: ind1, ind2, ind3, imaxrho 
    718775      !!---------------------------------------------------------------------- 
    719776!RB_CMEMS : declare here init for top       
     
    723780      ind2 = nn_hls + 2 + nbghostcells_x 
    724781      ind3 = nn_hls + 2 + nbghostcells_y_s 
     782      imaxrho = MAX(Agrif_irhox(), Agrif_irhoy()) 
    725783 
    726784      CALL agrif_declare_variable((/2,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/jpi,jpj,jpk,jptra+1/),trn_id) 
     
    735793      !----------------------------- 
    736794      CALL Agrif_Set_bc(trn_id,(/0,ind1-1/)) 
    737       CALL Agrif_Set_bc(trn_sponge_id,(/-nn_sponge_len*Agrif_irhox()-1,0/)) 
     795      CALL Agrif_Set_bc(trn_sponge_id,(/-nn_sponge_len*imaxrho-1,0/)) 
    738796 
    739797      ! 4. Update type 
  • NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/vremap.F90

    r12377 r13937  
    1 #undef PPR_LIB      /* USE PPR library */ 
     1#define PPR_LIB      /* USE PPR library */ 
    22MODULE vremap 
    33!$AGRIF_DO_NOT_TREAT 
     
    320320      ! specify methods 
    321321!      opts%edge_meth = p1e_method     ! 1st-order edge interp. 
     322!      opts%cell_meth = pcm_method 
    322323!      opts%cell_meth = plm_method     ! PLM method in cells 
    323324      opts%edge_meth = p3e_method     ! 3rd-order edge interp. 
     
    328329      ! limiter 
    329330!      opts%cell_lims = null_limit     ! no lim. 
     331!      opts%cell_lims = weno_limit 
    330332      opts%cell_lims = mono_limit     ! monotone limiter    
    331333  
     
    376378      DO jkout = 1, kjpk_out !  Loop over destination grid 
    377379         ! 
    378          IF     ( pzout(jkout) <=  pzin(  1    ) ) THEN  ; ptout(jkout,1:kn_var) = ptin(    1    ,1:kn_var)          
    379          ELSEIF ( pzout(jkout) >= pzin(kjpk_in) ) THEN   ; ptout(jkout,1:kn_var) = ptin( kjpk_in ,1:kn_var) 
     380         IF     ( pzout(jkout) <=  pzin(  1    ) ) THEN ! Surface extrapolation   
     381            DO jn = 1, kn_var  
     382               ptout(jkout,jn) = ptin(1 ,jn) + & 
     383                               & (pzout(jkout) - pzin(1)) / (pzin(2)    - pzin(1)) & 
     384                               &                          * (ptin(2,jn) - ptin(1,jn)) 
     385            END DO 
     386         ELSEIF ( pzout(jkout) >= pzin(kjpk_in) ) THEN ! Bottom extrapolation  
     387            DO jn = 1, kn_var  
     388               ptout(jkout,jn) = ptin(kjpk_in ,jn) + & 
     389                               & (pzout(jkout) - pzin(kjpk_in)) / (pzin(kjpk_in)    - pzin(kjpk_in-1)) & 
     390                               &                                * (ptin(kjpk_in,jn) - ptin(kjpk_in-1,jn)) 
     391            END DO 
    380392         ELSEIF ( ( pzout(jkout) > pzin(1) ).AND.( pzout(jkout) < pzin(kjpk_in) )) THEN 
    381393            DO jkin = 1, kjpk_in - 1 !  Loop over source grid 
Note: See TracChangeset for help on using the changeset viewer.