Skip to content
Snippets Groups Projects
agrif_oce_interp.F90 83 KiB
Newer Older
Guillaume Samson's avatar
Guillaume Samson committed
                        DO jk=2,N_in
                           z_in(jk) = z_in(jk-1) + 0.5_wp * ( h_in(jk) + h_in(jk-1) )
Guillaume Samson's avatar
Guillaume Samson committed
                        END DO
                        z_in(1:N_in) = z_in(1:N_in)  - ptab(ji,jj,k2,n2)
Guillaume Samson's avatar
Guillaume Samson committed
                     ENDIF                              

                     ! Output (Child) grid:
                     DO jk=1,N_out
                        h_out(jk) = e3t(ji,jj,jk,Krhs_a)
                     END DO
                     z_out(1) = 0.5_wp * e3w(ji,jj,1,Krhs_a) 
Guillaume Samson's avatar
Guillaume Samson committed
                     DO jk=2,N_out
                        z_out(jk) = z_out(jk-1) + e3w(ji,jj,jk,Krhs_a) 
Guillaume Samson's avatar
Guillaume Samson committed
                     END DO
                     IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out)  - ssh(ji,jj,Krhs_a)

                     IF( l_ini_child ) THEN
                        CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),          &
Guillaume Samson's avatar
Guillaume Samson committed
                                      &   z_out(1:N_out),N_in,N_out,jpts)  
                     ELSE 
                        CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),   &
Guillaume Samson's avatar
Guillaume Samson committed
                                      &   h_out(1:N_out),N_in,N_out,jpts)  
                     ENDIF
                  ENDIF
               END DO
            END DO
            Krhs_a = item
 
         ELSE
         
            IF ( Agrif_Parent(ln_zps) ) THEN ! Account for partial cells 
                                             ! linear vertical interpolation
               DO jj=j1,j2
                  DO ji=i1,i2
                     !
                     N_in  = mbkt(ji,jj)
                     N_out = mbkt(ji,jj)
                     z_in(1) = ptab(ji,jj,1,n2)
                     tabin(1,1:jpts) = ptab(ji,jj,1,1:jpts)
                     DO jk=2, N_in
                        z_in(jk) = ptab(ji,jj,jk,n2)
                        tabin(jk,1:jpts) = ptab(ji,jj,jk,1:jpts)
                     END DO
                     IF (.NOT.ln_linssh) z_in(1:N_in) = z_in(1:N_in) - ptab(ji,jj,k2,n2)
                     z_out(1) = 0.5_wp * e3w(ji,jj,1,Krhs_a)
Guillaume Samson's avatar
Guillaume Samson committed
                     DO jk=2, N_out
                        z_out(jk) = z_out(jk-1) + e3w(ji,jj,jk,Krhs_a) 
Guillaume Samson's avatar
Guillaume Samson committed
                     END DO
                     IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out) - ssh(ji,jj,Krhs_a)
                     CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),ptab(ji,jj,1:N_out,1:jpts), &
                                   &   z_out(1:N_out),N_in,N_out,jpts)  
                  END DO
               END DO
            ENDIF

            DO jn =1, jpts
               ts(i1:i2,j1:j2,1:jpkm1,jn,Krhs_a) = ptab(i1:i2,j1:j2,1:jpkm1,jn)*tmask(i1:i2,j1:j2,1:jpkm1)
Guillaume Samson's avatar
Guillaume Samson committed
            END DO
         ENDIF

      ENDIF
      !
   END SUBROUTINE interptsn

   
   SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE interpsshn  ***
      !!----------------------------------------------------------------------  
      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
      LOGICAL                         , INTENT(in   ) ::   before
      !
      !!----------------------------------------------------------------------  
      !
      IF( before) THEN
Guillaume Samson's avatar
Guillaume Samson committed
#if defined key_si3
         IF (l_ini_child.AND.(.NOT.(ln_rstart .OR. nn_iceini_file == 2))) THEN
            IF( ln_ice_embd ) THEN  
               ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Kmm_a)  &
                    &              + snwice_mass(i1:i2,j1:j2) * r1_rho0
            ELSE
               ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Kmm_a)  &
                    &              + rsshadj * tmask(i1:i2,j1:j2,1)
            ENDIF
         ELSE
            ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Kmm_a)
         ENDIF
#else
Guillaume Samson's avatar
Guillaume Samson committed
         ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Kmm_a)
Guillaume Samson's avatar
Guillaume Samson committed
#endif
Guillaume Samson's avatar
Guillaume Samson committed
      ELSE
         IF( l_ini_child ) THEN
            ssh(i1:i2,j1:j2,Krhs_a) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
         ELSE
            hbdy(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
         ENDIF
      ENDIF
      !
   END SUBROUTINE interpsshn

   
   SUBROUTINE interpsshn_frc( ptab, i1, i2, j1, j2, before )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE interpsshn  ***
      !!----------------------------------------------------------------------  
      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
      LOGICAL                         , INTENT(in   ) ::   before
      !
      !!----------------------------------------------------------------------  
      !
      IF( before) THEN
         ptab(i1:i2,j1:j2) = ssh_frc(i1:i2,j1:j2)
      ELSE
         ssh_frc(i1:i2,j1:j2) = ptab(i1:i2,j1:j2)
      ENDIF
      !
   END SUBROUTINE interpsshn_frc


   SUBROUTINE interp_tmask_agrif( ptab, i1, i2, j1, j2, before )
      !!----------------------------------------------------------------------
      !!               ***  ROUTINE interp_tmask_agrif  ***
      !!
      !!               set tmask_agrif = 0 over ghost points 
      !!
      !!----------------------------------------------------------------------  
      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
      LOGICAL                         , INTENT(in   ) ::   before
      !
      !!----------------------------------------------------------------------  
      !
      IF(.NOT.before) THEN
         tmask_agrif(i1:i2,j1:j2) = 0._wp 
      ENDIF
      !
   END SUBROUTINE interp_tmask_agrif


Guillaume Samson's avatar
Guillaume Samson committed
   SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
      !!----------------------------------------------------------------------
      !!                  *** ROUTINE interpun ***
      !!---------------------------------------------    
      !!
      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
      LOGICAL, INTENT(in) :: before
      !!
      INTEGER :: ji,jj,jk
Guillaume Samson's avatar
Guillaume Samson committed
      ! vertical interpolation:
      REAL(wp), DIMENSION(i1:i2,j1:j2) :: zsshu
Guillaume Samson's avatar
Guillaume Samson committed
      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in, z_in
      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out
Guillaume Samson's avatar
Guillaume Samson committed
      !!---------------------------------------------    
      !
      IF (before) THEN 

         item = Kmm_a
         IF( l_ini_child )   Kmm_a = Kbb_a     

Guillaume Samson's avatar
Guillaume Samson committed
            DO jj=j1,j2
               DO ji=i1,i2
                  ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a)*umask(ji,jj,jk)) 
                  !!IF( l_vremap .OR. l_ini_child) THEN
                  !!   ! Interpolate thicknesses (masked for subsequent extrapolation)
                  !!   ptab(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
                  !!ENDIF
Guillaume Samson's avatar
Guillaume Samson committed
               END DO
            END DO
         END DO

         Kmm_a = item
Guillaume Samson's avatar
Guillaume Samson committed
      ELSE
         zrhoy = Agrif_rhoy()

Guillaume Samson's avatar
Guillaume Samson committed

            zsshu(i1:i2,j1:j2) = 0._wp 
            
            IF ( .NOT.ln_linssh ) THEN
               zsshu(i1:i2,j1:j2) = hu(i1:i2,j1:j2,Krhs_a) - hu_0(i1:i2,j1:j2)   
            ENDIF   
Guillaume Samson's avatar
Guillaume Samson committed

            DO ji=i1,i2
               DO jj=j1,j2
                  uu(ji,jj,:,Krhs_a) = 0._wp
Guillaume Samson's avatar
Guillaume Samson committed
                  N_out = mbku(ji,jj)
                  IF (N_in*N_out > 0) THEN
Guillaume Samson's avatar
Guillaume Samson committed
                     DO jk=1,N_in
                        h_in(jk)  = e3u0_parent(ji,jj,jk) * & 
                             &       (1._wp + zsshu(ji,jj)/(hu0_parent(ji,jj)*ssumask(ji,jj) + 1._wp - ssumask(ji,jj)))
                        tabin(jk) = ptab(ji,jj,jk,1) / (e2u(ji,jj)*zrhoy*h_in(jk))
                     END DO
Guillaume Samson's avatar
Guillaume Samson committed
                     
                     DO jk=1, N_out
                        h_out(jk) = e3u(ji,jj,jk,Krhs_a)
                     END DO

                     IF( l_ini_child ) THEN
                        z_in(1) = 0.5_wp * h_in(1)
                        DO jk=2,N_in
                           z_in(jk) = z_in(jk-1) + 0.5_wp * (h_in(jk)+h_in(jk-1))
                        END DO
                        !
                        z_out(1) = 0.5_wp * h_out(1)
                        DO jk=2,N_out
                           z_out(jk) = z_out(jk-1) + 0.5_wp * (h_out(jk-1) + h_out(jk)) 
                        END DO  

                        CALL remap_linear       (tabin(1:N_in),z_in(1:N_in),uu(ji,jj,1:N_out,Krhs_a),z_out(1:N_out),N_in,N_out,1)
                     ELSE
                        CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),uu(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out,1)
                     ENDIF   
                  ENDIF
Guillaume Samson's avatar
Guillaume Samson committed
               END DO
            END DO
         ELSE
            DO jk = 1, jpkm1
               uu(i1:i2,j1:j2,jk,Krhs_a) = ptab(i1:i2,j1:j2,jk,1) / ( zrhoy * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Krhs_a) )
Guillaume Samson's avatar
Guillaume Samson committed
            END DO
         ENDIF

      ENDIF
      ! 
   END SUBROUTINE interpun

   
   SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
      !!----------------------------------------------------------------------
      !!                  *** ROUTINE interpvn ***
      !!----------------------------------------------------------------------
      !
      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
      LOGICAL, INTENT(in) :: before
      !
      INTEGER :: ji,jj,jk
      REAL(wp) :: zrhox
      ! vertical interpolation:
      REAL(wp), DIMENSION(i1:i2,j1:j2) :: zsshv
Guillaume Samson's avatar
Guillaume Samson committed
      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in, z_in
      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out
      INTEGER  :: N_in, N_out, item
      !!---------------------------------------------    
      !      
      IF (before) THEN   

         item = Kmm_a
         IF( l_ini_child )   Kmm_a = Kbb_a     
       
Guillaume Samson's avatar
Guillaume Samson committed
            DO jj=j1,j2
               DO ji=i1,i2
                  ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)*vmask(ji,jj,jk))
                  !!IF( l_vremap .OR. l_ini_child) THEN
                  !!   ! Interpolate thicknesses (masked for subsequent extrapolation)
                  !!   ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
                  !!ENDIF
Guillaume Samson's avatar
Guillaume Samson committed
               END DO
            END DO
         END DO

Guillaume Samson's avatar
Guillaume Samson committed

      ELSE       
         zrhox = Agrif_rhox()

Guillaume Samson's avatar
Guillaume Samson committed

            zsshv(i1:i2,j1:j2) = 0._wp 
            
            IF ( .NOT.ln_linssh ) THEN
               zsshv(i1:i2,j1:j2) = hv(i1:i2,j1:j2,Krhs_a) - hv_0(i1:i2,j1:j2)   
            ENDIF   
Guillaume Samson's avatar
Guillaume Samson committed

Guillaume Samson's avatar
Guillaume Samson committed
                  vv(ji,jj,:,Krhs_a) = 0._wp
Guillaume Samson's avatar
Guillaume Samson committed
                  N_out = mbkv(ji,jj)
                  IF (N_in*N_out > 0) THEN

                     DO jk=1,N_in
                        h_in(jk)  = e3v0_parent(ji,jj,jk) * & 
                             &       (1._wp + zsshv(ji,jj)/(hv0_parent(ji,jj)*ssvmask(ji,jj) + 1._wp - ssvmask(ji,jj)))
                        tabin(jk) = ptab(ji,jj,jk,1) / (e1v(ji,jj)*zrhox*h_in(jk))
Guillaume Samson's avatar
Guillaume Samson committed
                     END DO
Guillaume Samson's avatar
Guillaume Samson committed
                        h_out(jk) = e3v(ji,jj,jk,Krhs_a)
                     END DO

                     IF( l_ini_child ) THEN
                        z_in(1) = 0.5_wp * h_in(1)
                        DO jk=2,N_in
                           z_in(jk) = z_in(jk-1) + 0.5_wp * (h_in(jk)+h_in(jk-1))
                        END DO
                        !
                        z_out(1) = 0.5_wp * h_out(1)
                        DO jk=2,N_out
                           z_out(jk) = z_out(jk-1) + 0.5_wp * (h_out(jk-1) + h_out(jk)) 
                        END DO  

Guillaume Samson's avatar
Guillaume Samson committed
                        CALL remap_linear       (tabin(1:N_in),z_in(1:N_in),vv(ji,jj,1:N_out,Krhs_a),z_out(1:N_out),N_in,N_out,1)
                     ELSE
                        CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),vv(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out,1)
                     ENDIF   
                  ENDIF
               END DO
            END DO
         ELSE
            DO jk = 1, jpkm1
               vv(i1:i2,j1:j2,jk,Krhs_a) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v(i1:i2,j1:j2,jk,Krhs_a) )
            END DO
         ENDIF
      ENDIF
      !        
   END SUBROUTINE interpvn

   SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before)
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE interpunb  ***
      !!----------------------------------------------------------------------  
      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
      LOGICAL                         , INTENT(in   ) ::   before
      !
      INTEGER  ::   ji, jj
      REAL(wp) ::   zrhoy, zrhot, zt0, zt1, ztcoeff
      !!----------------------------------------------------------------------  
      !
      IF( before ) THEN 
         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu(i1:i2,j1:j2,Kmm_a) * uu_b(i1:i2,j1:j2,Kmm_a)
      ELSE
         zrhoy = Agrif_Rhoy()
         zrhot = Agrif_rhot()
         ! Time indexes bounds for integration
         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot      
         ! 
         DO ji = i1, i2
            DO jj = j1, j2
               IF ( utint_stage(ji,jj) < (bdy_tinterp + 1) ) THEN
                  IF    ( utint_stage(ji,jj) == 1  ) THEN
                     ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
                        &               - zt0**2._wp * (       zt0 - 1._wp)        )
                  ELSEIF( utint_stage(ji,jj) == 2  ) THEN
                     ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
                        &               - zt0        * (       zt0 - 1._wp)**2._wp )
                  ELSEIF( utint_stage(ji,jj) == 0  ) THEN                
                     ztcoeff = 1._wp
                  ELSE
                     ztcoeff = 0._wp
                  ENDIF
                  !   
                  ubdy(ji,jj) = ubdy(ji,jj) + ztcoeff * ptab(ji,jj)
                  !            
                  IF (( utint_stage(ji,jj) == 2 ).OR.( utint_stage(ji,jj) == 0 )) THEN
                     ubdy(ji,jj) = ubdy(ji,jj) / (zrhoy*e2u(ji,jj)) * umask(ji,jj,1)
                  ENDIF
                  !
                  utint_stage(ji,jj) = utint_stage(ji,jj) + 1
               ENDIF
            END DO
         END DO
      END IF
      ! 
   END SUBROUTINE interpunb


   SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE interpvnb  ***
      !!----------------------------------------------------------------------  
      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
      LOGICAL                         , INTENT(in   ) ::   before
      !
      INTEGER  ::   ji, jj
      REAL(wp) ::   zrhox, zrhot, zt0, zt1, ztcoeff   
      !!----------------------------------------------------------------------  
      ! 
      IF( before ) THEN 
         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv(i1:i2,j1:j2,Kmm_a) * vv_b(i1:i2,j1:j2,Kmm_a)
      ELSE
         zrhox = Agrif_Rhox()
         zrhot = Agrif_rhot()
         ! Time indexes bounds for integration
         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot 
         !     
         DO ji = i1, i2
            DO jj = j1, j2
               IF ( vtint_stage(ji,jj) < (bdy_tinterp + 1) ) THEN
                  IF    ( vtint_stage(ji,jj) == 1  ) THEN
                     ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
                        &               - zt0**2._wp * (       zt0 - 1._wp)        )
                  ELSEIF( vtint_stage(ji,jj) == 2  ) THEN
                     ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
                        &               - zt0        * (       zt0 - 1._wp)**2._wp )
                  ELSEIF( vtint_stage(ji,jj) == 0  ) THEN                
                     ztcoeff = 1._wp
                  ELSE
                     ztcoeff = 0._wp
                  ENDIF
                  !   
                  vbdy(ji,jj) = vbdy(ji,jj) + ztcoeff * ptab(ji,jj)
                  !            
                  IF (( vtint_stage(ji,jj) == 2 ).OR.( vtint_stage(ji,jj) == 0 )) THEN
                     vbdy(ji,jj) = vbdy(ji,jj) / (zrhox*e1v(ji,jj)) * vmask(ji,jj,1)
                  ENDIF
                  !
                  vtint_stage(ji,jj) = vtint_stage(ji,jj) + 1
               ENDIF
            END DO
         END DO          
      ENDIF
      !
   END SUBROUTINE interpvnb


   SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE interpub2b  ***
      !!----------------------------------------------------------------------  
      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
      LOGICAL                         , INTENT(in   ) ::   before
      !
      INTEGER  ::   ji,jj
      REAL(wp) ::   zrhot, zt0, zt1, zat
      !!----------------------------------------------------------------------  
      IF( before ) THEN
!         IF ( ln_bt_fw ) THEN
            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2)
!         ELSE
!            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * un_adv(i1:i2,j1:j2)
!         ENDIF
      ELSE
         zrhot = Agrif_rhot()
         ! Time indexes bounds for integration
         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
         ! Polynomial interpolation coefficients:
         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
         !
         ubdy(i1:i2,j1:j2) = zat * ptab(i1:i2,j1:j2) 
         !
         ! Update interpolation stage:
         utint_stage(i1:i2,j1:j2) = 1
      ENDIF
      ! 
   END SUBROUTINE interpub2b
   
   SUBROUTINE interpub2b_const( ptab, i1, i2, j1, j2, before )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE interpub2b_const  ***
      !!----------------------------------------------------------------------  
      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
      LOGICAL                         , INTENT(in   ) ::   before
      !
      REAL(wp) :: zrhoy
      !!----------------------------------------------------------------------  
      IF( before ) THEN
!         IF ( ln_bt_fw ) THEN
            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2) &
                                * umask(i1:i2,j1:j2,1)
Guillaume Samson's avatar
Guillaume Samson committed
!         ELSE
!            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * un_adv(i1:i2,j1:j2)
!         ENDIF
      ELSE
         zrhoy = Agrif_Rhoy()
         !
         ubdy(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) & 
                           & / (zrhoy*e2u(i1:i2,j1:j2)) * umask(i1:i2,j1:j2,1)
         !
      ENDIF
      ! 
   END SUBROUTINE interpub2b_const


   SUBROUTINE ub2b_cor( ptab, i1, i2, j1, j2, before )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE ub2b_cor  ***
      !!----------------------------------------------------------------------  
      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
      LOGICAL                         , INTENT(in   ) ::   before
      !
      INTEGER  :: ji, jj
      INTEGER  :: imin, imax, jmin, jmax
      REAL(wp) :: zrhox, zrhoy, zx
      !!----------------------------------------------------------------------  
      IF( before ) THEN
         ptab(:,:) = 0._wp
         imin = MAX(i1, 2) ; imax = MIN(i2, jpi-1)
         jmin = MAX(j1, 2) ; jmax = MIN(j2, jpj-1)
         DO ji=imin,imax
            DO jj=jmin,jmax
               ptab(ji,jj) = 0.25_wp *(vmask(ji,jj  ,1)                        & 
                           &       * ( vb2_b(ji+1,jj  )*e1v(ji+1,jj  )         & 
                           &          -vb2_b(ji-1,jj  )*e1v(ji-1,jj  ) )       &
                           &          -vmask(ji,jj-1,1)                        & 
                           &       * ( vb2_b(ji+1,jj-1)*e1v(ji+1,jj-1)         &
                           &          -vb2_b(ji-1,jj-1)*e1v(ji-1,jj-1) ) )      
Guillaume Samson's avatar
Guillaume Samson committed
            END DO
         END DO 
      ELSE
         !
         zrhox = Agrif_Rhox() 
         zrhoy = Agrif_Rhoy()
         DO ji=i1,i2
            DO jj=j1,j2
               IF (utint_stage(ji,jj)==0) THEN 
                  zx = 2._wp*MOD(ABS(mig0(ji)-nbghostcells_x_w), INT(Agrif_Rhox()))/zrhox - 1._wp  
                  ubdy(ji,jj) = ubdy(ji,jj) + 0.25_wp*(1._wp-zx*zx) * ptab(ji,jj) & 
                              &         / zrhoy *r1_e2u(ji,jj) * umask(ji,jj,1) 
                  utint_stage(ji,jj) = 1 
               ENDIF
            END DO
         END DO 
         !
      ENDIF
      ! 
   END SUBROUTINE ub2b_cor


   SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE interpvb2b  ***
      !!----------------------------------------------------------------------  
      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
      LOGICAL                         , INTENT(in   ) ::   before
      !
      INTEGER ::   ji,jj
      REAL(wp) ::   zrhot, zt0, zt1, zat
      !!----------------------------------------------------------------------  
      !
      IF( before ) THEN
!         IF ( ln_bt_fw ) THEN
            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2)
!         ELSE
!            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2)
!         ENDIF
      ELSE      
         zrhot = Agrif_rhot()
         ! Time indexes bounds for integration
         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
         ! Polynomial interpolation coefficients:
         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
         !
         vbdy(i1:i2,j1:j2) = zat * ptab(i1:i2,j1:j2)
         !
         ! update interpolation stage:
         vtint_stage(i1:i2,j1:j2) = 1
      ENDIF
      !      
   END SUBROUTINE interpvb2b


   SUBROUTINE interpvb2b_const( ptab, i1, i2, j1, j2, before )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE interpub2b_const  ***
      !!----------------------------------------------------------------------  
      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
      LOGICAL                         , INTENT(in   ) ::   before
      !
      REAL(wp) :: zrhox
      !!----------------------------------------------------------------------  
      IF( before ) THEN
!         IF ( ln_bt_fw ) THEN
            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2) &
                                * vmask(i1:i2,j1:j2,1)
Guillaume Samson's avatar
Guillaume Samson committed
!         ELSE
!            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2)
!         ENDIF
      ELSE
         zrhox = Agrif_Rhox()
         !
         vbdy(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) &
                           & / (zrhox*e1v(i1:i2,j1:j2)) * vmask(i1:i2,j1:j2,1)
         !
      ENDIF
      ! 
   END SUBROUTINE interpvb2b_const

 
   SUBROUTINE vb2b_cor( ptab, i1, i2, j1, j2, before )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE vb2b_cor  ***
      !!----------------------------------------------------------------------  
      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
      LOGICAL                         , INTENT(in   ) ::   before
      !
      INTEGER  :: ji, jj
      INTEGER  :: imin, imax, jmin, jmax
      REAL(wp) :: zrhox, zrhoy, zy, zslope1, zslope2
Guillaume Samson's avatar
Guillaume Samson committed
      !!----------------------------------------------------------------------  
      IF( before ) THEN
         ptab(:,:) = 0._wp
         imin = MAX(i1, 2) ; imax = MIN(i2, jpi-1)
         jmin = MAX(j1, 2) ; jmax = MIN(j2, jpj-1)
         DO ji=imin,imax
            DO jj=jmin,jmax
               ptab(ji,jj) = 0.25_wp *(umask(ji  ,jj,1)                      & 
                           &       * ( ub2_b(ji  ,jj+1)*e2u(ji  ,jj+1)       & 
                           &          -ub2_b(ji  ,jj-1)*e2u(ji  ,jj-1) )     &
                           &          -umask(ji-1,jj,1)                      & 
                           &       * ( ub2_b(ji-1,jj+1)*e2u(ji-1,jj+1)       &
                           &          -ub2_b(ji-1,jj-1)*e2u(ji-1,jj-1) ) )   
Guillaume Samson's avatar
Guillaume Samson committed
            END DO
         END DO 
      ELSE
         !
         zrhox = Agrif_Rhox() 
         zrhoy = Agrif_Rhoy()
         DO ji=i1,i2
            DO jj=j1,j2
               IF (vtint_stage(ji,jj)==0) THEN 
                  zy = 2._wp*MOD(ABS(mjg0(jj)-nbghostcells_y_s), INT(Agrif_Rhoy()))/zrhoy - 1._wp  
                  vbdy(ji,jj) = vbdy(ji,jj) + 0.25_wp*(1._wp-zy*zy) * ptab(ji,jj) & 
                              &         / zrhox * r1_e1v(ji,jj) * vmask(ji,jj,1) 
                  vtint_stage(ji,jj) = 1 
               ENDIF
            END DO
         END DO 
         !
      ENDIF
      ! 
   END SUBROUTINE vb2b_cor


   SUBROUTINE interpe3t0_vremap( ptab, i1, i2, j1, j2, k1, k2, before )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE interpe3t0_vremap  ***
      !!----------------------------------------------------------------------  
      INTEGER                              , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
      LOGICAL                              , INTENT(in   ) :: before
      !
      INTEGER :: ji, jj, jk
      REAL(wp) :: zh
      !!----------------------------------------------------------------------  
      !    
      IF( before ) THEN
         IF ( ln_zps ) THEN
            DO jk = k1, k2
               DO jj = j1, j2
                  DO ji = i1, i2
                     ptab(ji, jj, jk) = e3t_1d(jk)
                  END DO
               END DO
            END DO
         ELSE
            ptab(i1:i2,j1:j2,k1:k2) = e3t_0(i1:i2,j1:j2,k1:k2)
         ENDIF
      ELSE
         !
         DO jk = k1, k2
            DO jj = j1, j2
               DO ji = i1, i2
                  e3t0_parent(ji,jj,jk) = ptab(ji,jj,jk)
               END DO
            END DO
         END DO

         ! Retrieve correct scale factor at the bottom:
         DO jj = j1, j2
            DO ji = i1, i2
               IF ( mbkt_parent(ji,jj) > 1 ) THEN
                  zh = 0._wp
                  DO jk = 1, mbkt_parent(ji, jj)-1
                     zh = zh + e3t0_parent(ji,jj,jk)
                  END DO
                  e3t0_parent(ji,jj,mbkt_parent(ji,jj)) = ht0_parent(ji, jj) - zh
               ENDIF  
            END DO
         END DO
         
      ENDIF
      ! 
   END SUBROUTINE interpe3t0_vremap


   SUBROUTINE interpglamt( ptab, i1, i2, j1, j2, before )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE interpglamt  ***
      !!----------------------------------------------------------------------  
      INTEGER                        , INTENT(in   ) :: i1, i2, j1, j2
      REAL(wp),DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
      LOGICAL                        , INTENT(in   ) :: before
      !
      INTEGER :: ji, jj, jk
      REAL(wp):: ztst
      !!----------------------------------------------------------------------  
      !    
      IF( before ) THEN
         ptab(i1:i2,j1:j2) = glamt(i1:i2,j1:j2)
      ELSE
         ztst = MAXVAL(ABS(glamt(i1:i2,j1:j2)))*1.e-4
         DO jj = j1, j2
            DO ji = i1, i2
               IF( ABS( ptab(ji,jj) - glamt(ji,jj) ) > ztst ) THEN
                  WRITE(numout,*) ' Agrif error for glamt: parent, child, i, j ', ptab(ji,jj), glamt(ji,jj), mig0(ji), mig0(jj)
!                  kindic_agr = kindic_agr + 1
               ENDIF
            END DO
         END DO
      ENDIF
      ! 
   END SUBROUTINE interpglamt


   SUBROUTINE interpgphit( ptab, i1, i2, j1, j2, before )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE interpgphit  ***
      !!----------------------------------------------------------------------  
      INTEGER                        , INTENT(in   ) :: i1, i2, j1, j2
      REAL(wp),DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
      LOGICAL                        , INTENT(in   ) :: before
      !
      INTEGER :: ji, jj, jk
      REAL(wp):: ztst
      !!----------------------------------------------------------------------  
      !    
      IF( before ) THEN
         ptab(i1:i2,j1:j2) = gphit(i1:i2,j1:j2)
      ELSE
         ztst = MAXVAL(ABS(gphit(i1:i2,j1:j2)))*1.e-4
         DO jj = j1, j2
            DO ji = i1, i2
               IF( ABS( ptab(ji,jj) - gphit(ji,jj) ) > ztst ) THEN
                  WRITE(numout,*) ' Agrif error for gphit: parent, child, i, j ', ptab(ji,jj), gphit(ji,jj), mig0(ji), mig0(jj)
!                  kindic_agr = kindic_agr + 1
               ENDIF
            END DO
         END DO
      ENDIF
      ! 
   END SUBROUTINE interpgphit


   SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE interavm  ***
      !!----------------------------------------------------------------------  
      INTEGER                                    , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, m1, m2
      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) ::   ptab
      LOGICAL                                    , INTENT(in   ) ::   before
      !
      INTEGER  :: ji, jj, jk
      INTEGER  :: N_in, N_out
      REAL(wp), DIMENSION(k1:k2) :: tabin, z_in
      REAL(wp), DIMENSION(1:jpk) :: z_out
      !!----------------------------------------------------------------------  
      !      
      IF (before) THEN         
         DO jk=k1,k2
            DO jj=j1,j2
              DO ji=i1,i2
                    ptab(ji,jj,jk,1) = avm_k(ji,jj,jk)
              END DO
           END DO
         END DO

         IF( l_vremap ) THEN
            ! Interpolate interfaces 
            ! Warning: these are masked, hence extrapolated prior interpolation.
            DO jk=k1,k2
               DO jj=j1,j2
                  DO ji=i1,i2
                      ptab(ji,jj,jk,2) = tmask(ji,jj,jk) * gdepw(ji,jj,jk,Kmm_a)
                  END DO
               END DO
            END DO
        
           ! Save ssh at last level:
            IF (.NOT.ln_linssh) THEN
               ptab(i1:i2,j1:j2,k2,2) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 
            ELSE
               ptab(i1:i2,j1:j2,k2,2) = 0._wp
            END IF      
          ENDIF

      ELSE 

         IF( l_vremap ) THEN
            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
            avm_k(i1:i2,j1:j2,1:jpkm1) = 0._wp
               
            DO jj = j1, j2
               DO ji =i1, i2
                  N_in = mbkt_parent(ji,jj)
                  N_out = mbkt(ji,jj)
                  IF (N_in*N_out > 0) THEN
                     DO jk = 1, N_in  ! Parent vertical grid               
                        z_in(jk)  = ptab(ji,jj,jk,2) - ptab(ji,jj,k2,2)
                        tabin(jk) = ptab(ji,jj,jk,1)
                     END DO
                     DO jk = 1, N_out        ! Child vertical grid
                        z_out(jk) = gdepw(ji,jj,jk,Kmm_a) - ssh(ji,jj,Kmm_a)
                     END DO
                     IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out)  - ssh(ji,jj,Kmm_a)

                     CALL remap_linear(tabin(1:N_in),z_in(1:N_in),avm_k(ji,jj,1:N_out),z_out(1:N_out),N_in,N_out,1)
                  ENDIF
               END DO
            END DO
         ELSE
            avm_k(i1:i2,j1:j2,1:jpkm1) = ptab (i1:i2,j1:j2,1:jpkm1,1)
         ENDIF
      ENDIF
      !
   END SUBROUTINE interpavm

   
   SUBROUTINE interpmbkt( ptab, i1, i2, j1, j2, before )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE interpmbkt  ***
      !!----------------------------------------------------------------------  
      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
      LOGICAL                         , INTENT(in   ) ::   before
      !
      !!----------------------------------------------------------------------  
      !
      IF( before) THEN
         ptab(i1:i2,j1:j2) = REAL(mbkt(i1:i2,j1:j2),wp)
      ELSE
         mbkt_parent(i1:i2,j1:j2) = NINT(ptab(i1:i2,j1:j2))
      ENDIF
      !
   END SUBROUTINE interpmbkt

   
   SUBROUTINE interpht0( ptab, i1, i2, j1, j2, before )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE interpht0  ***
      !!----------------------------------------------------------------------  
      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
      LOGICAL                         , INTENT(in   ) ::   before
      !
      !!----------------------------------------------------------------------  
      !
      IF( before) THEN
         ptab(i1:i2,j1:j2) = ht_0(i1:i2,j1:j2) * ssmask(i1:i2,j1:j2)
Guillaume Samson's avatar
Guillaume Samson committed
      ELSE
         ht0_parent(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) * ssmask(i1:i2,j1:j2)
Guillaume Samson's avatar
Guillaume Samson committed
      ENDIF
      !
   END SUBROUTINE interpht0


   SUBROUTINE interp_e1e2t_frac(tabres, i1, i2, j1, j2, before )
      !
      !!----------------------------------------------------------------------
      !!               *** ROUTINE interp_e1e2t_frac ***
      !!----------------------------------------------------------------------
      INTEGER                        , INTENT(in   ) :: i1, i2, j1, j2
      REAL(wp),DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
      LOGICAL                        , INTENT(in   ) :: before
      !!
      !!----------------------------------------------------------------------

      IF (before) THEN
         tabres(i1:i2,j1:j2) = e1e2t(i1:i2,j1:j2)
      ELSE
         WHERE (tabres(i1:i2,j1:j2)/=0._wp)
            e1e2t_frac(i1:i2,j1:j2) = e1e2t(i1:i2,j1:j2) &
                 & / tabres(i1:i2,j1:j2) * Agrif_Rhox() * Agrif_Rhoy()
         ELSEWHERE
            e1e2t_frac(i1:i2,j1:j2) = 1._wp
         END WHERE
      ENDIF
      !
   END SUBROUTINE interp_e1e2t_frac


   SUBROUTINE interp_e2u_frac(tabres, i1, i2, j1, j2, before )
      !
      !!----------------------------------------------------------------------
      !!               *** ROUTINE interp_e2u_frac ***
      !!----------------------------------------------------------------------
      INTEGER                        , INTENT(in   ) :: i1, i2, j1, j2
      REAL(wp),DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
      LOGICAL                        , INTENT(in   ) :: before
      !!
      !!----------------------------------------------------------------------

      IF (before) THEN
         tabres(i1:i2,j1:j2) = e2u(i1:i2,j1:j2)
      ELSE
         WHERE (tabres(i1:i2,j1:j2)/=0._wp)
            e2u_frac(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) &
                 & / tabres(i1:i2,j1:j2) * Agrif_Rhoy()
         ELSE WHERE
            e2u_frac(i1:i2,j1:j2) = 1._wp
         END WHERE
      ENDIF
      !
   END SUBROUTINE interp_e2u_frac


   SUBROUTINE interp_e1v_frac(tabres, i1, i2, j1, j2, before )
      !
      !!----------------------------------------------------------------------
      !!               *** ROUTINE interp_e1v_frac ***
      !!----------------------------------------------------------------------
      INTEGER                        , INTENT(in   ) :: i1, i2, j1, j2
      REAL(wp),DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
      LOGICAL                        , INTENT(in   ) :: before
      !!
      !!----------------------------------------------------------------------

      IF (before) THEN
         tabres(i1:i2,j1:j2) = e1v(i1:i2,j1:j2)
      ELSE
         WHERE (tabres(i1:i2,j1:j2)/=0._wp)
            e1v_frac(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) &
                 & / tabres(i1:i2,j1:j2) * Agrif_Rhox()
         ELSE WHERE
            e1v_frac(i1:i2,j1:j2) = 1._wp
         END WHERE
      ENDIF
      !
   END SUBROUTINE interp_e1v_frac


Guillaume Samson's avatar
Guillaume Samson committed
   SUBROUTINE Agrif_check_bat( iindic )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE Agrif_check_bat  ***
      !!----------------------------------------------------------------------  
      INTEGER, INTENT(inout) ::   iindic
      !!
      INTEGER :: ji, jj, jk
      INTEGER  :: istart, iend, jstart, jend, ispon
      !!----------------------------------------------------------------------  
      !
      !
      ! --- West --- !
      IF(lk_west) THEN
         ispon  = (nn_sponge_len+2) * Agrif_irhox()
Guillaume Samson's avatar
Guillaume Samson committed
         istart = nn_hls + 2                                  ! halo + land + 1
         iend   = nn_hls + nbghostcells + ispon           ! halo + land + nbghostcells + sponge
         jstart = nn_hls + 2 
         jend   = jpjglo - nn_hls - 1 
         DO ji = mi0(istart), mi1(iend)
            DO jj = mj0(jstart), mj1(jend)
               IF ( ABS(ht0_parent(ji,jj)-ht_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
               IF ( .NOT.ln_vert_remap) THEN
                  DO jk = 1, jpkm1
                     IF ( ABS(e3t0_parent(ji,jj,jk)-e3t_0(ji,jj,jk))*tmask(ji,jj,jk) > 1.e-3 ) iindic = iindic + 1
                  END DO 
               ENDIF
            END DO
            DO jj = mj0(jstart), mj1(jend-1)
               IF ( ABS(hv0_parent(ji,jj)-hv_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
               IF ( .NOT.ln_vert_remap) THEN
                  DO jk = 1, jpkm1
                     IF ( ABS(e3v0_parent(ji,jj,jk)-e3v_0(ji,jj,jk))*vmask(ji,jj,jk) > 1.e-3 ) iindic = iindic + 1
                  END DO 
               ENDIF
            END DO
         END DO
         DO ji = mi0(istart), mi1(iend-1)
            DO jj = mj0(jstart), mj1(jend)
               IF ( ABS(hu0_parent(ji,jj)-hu_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
               IF ( .NOT.ln_vert_remap) THEN
                  DO jk = 1, jpkm1
                     IF ( ABS(e3u0_parent(ji,jj,jk)-e3u_0(ji,jj,jk))*umask(ji,jj,jk) > 1.e-3 ) iindic = iindic + 1
                  END DO 
               ENDIF
            END DO
         END DO
      ENDIF
      !
      ! --- East --- !
      IF(lk_east) THEN
         ispon  = (nn_sponge_len+2) * Agrif_irhox() 
Guillaume Samson's avatar
Guillaume Samson committed
         istart = jpiglo - ( nn_hls + nbghostcells + ispon -1 )  ! halo + land + nbghostcells + sponge - 1
         iend   = jpiglo - nn_hls - 1                            ! halo + land + 1                     - 1
         jstart = nn_hls + 2 
         jend   = jpjglo - nn_hls - 1
         DO ji = mi0(istart), mi1(iend)