Skip to content
Snippets Groups Projects
crsdom.F90 112 KiB
Newer Older
Guillaume Samson's avatar
Guillaume Samson committed
MODULE crsdom
   !!===================================================================
   !!                  ***  crs.F90 ***
   !!  Purpose: Interface for calculating quantities from a  
   !!           higher-resolution grid for the coarse grid.
   !!
   !!  Method:  Given the user-defined reduction factor, 
   !!           the averaging bins are set: 
   !!           - nn_binref = 0, starting from the north 
   !!           to the south in the model interior domain, 
   !!           in this way the north fold and redundant halo cells  
   !!           could be handled in a consistent manner and 
   !!           the irregularities of bin size can be handled
   !!           more naturally by the presence of land
   !!           in the southern boundary.  Thus the southernmost bin 
   !!           could be of an irregular bin size.
   !!           Information on the parent grid is retained, specifically,
   !!           each coarse grid cell's volume and ocean surface 
   !!           at the faces, relative to the parent grid.
   !!           - nn_binref = 1 (not yet available), starting
   !!           at a centralized bin at the equator, being only
   !!           truly centered for odd-numbered j-direction reduction 
   !!           factors.
   !!  References:  Aumont, O., J.C. Orr, D. Jamous, P. Monfray
   !!               O. Marti and G. Madec, 1998. A degradation 
   !!               approach to accelerate simulations to steady-state 
   !!               in a 3-D tracer transport model of the global ocean.
   !!               Climate Dynamics, 14:101-116.
   !!  History: 
   !!       Original.   May 2012.  (J. Simeon, C. Calone, G. Madec, C. Ethe)
   !!===================================================================
   USE dom_oce        ! ocean space and time domain
   USE crs            ! domain for coarse grid
   !
   USE in_out_manager 
   USE par_kind
   USE crslbclnk
   USE lib_mpp

   IMPLICIT NONE

   PRIVATE 

   PUBLIC crs_dom_ope
   PUBLIC crs_dom_e3, crs_dom_sfc, crs_dom_msk, crs_dom_hgr, crs_dom_coordinates
   PUBLIC crs_dom_facvol, crs_dom_def, crs_dom_bat

   INTERFACE crs_dom_ope
      MODULE PROCEDURE crs_dom_ope_3d, crs_dom_ope_2d
   END INTERFACE

   REAL(wp) :: r_inf = 1e+36

   !!----------------------------------------------------------------------
   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
   !! $Id: crsdom.F90 14433 2021-02-11 08:06:49Z smasson $
   !! Software governed by the CeCILL license (see ./LICENSE)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE crs_dom_msk
      
      INTEGER  ::  ji, jj, jk                   ! dummy loop indices
      INTEGER  ::  ijie,ijis,ijje,ijjs,ij,je_2
      REAL(wp) ::  zmask
      
      ! Initialize

      tmask_crs(:,:,:) = 0.0
      vmask_crs(:,:,:) = 0.0
      umask_crs(:,:,:) = 0.0
      fmask_crs(:,:,:) = 0.0
  
            
      IF( Njs0_crs == 1 .AND. ( ( mje_crs(2) - mjs_crs(2) ) < 2 ) ) THEN     !!cc bande du sud style ORCA2
         IF( mje_crs(2) - mjs_crs(2) == 1 ) THEN
            je_2 = mje_crs(2)   ;  ij = je_2
         ENDIF
      ELSE
         je_2 = mje_crs(2)      ;  ij = mjs_crs(2) 
      ENDIF
      DO jk = 1, jpkm1
         DO ji = 2, Nie0_crs  
            ijis = mis_crs(ji)  ;  ijie = mie_crs(ji)    
            !          
            zmask = 0.0
            zmask = SUM( tmask(ijis:ijie,ij:je_2,jk) ) 
            IF ( zmask > 0.0 ) tmask_crs(ji,2,jk) = 1.0_wp
               
            zmask = 0.0
            zmask = SUM( vmask(ijis:ijie,je_2     ,jk) )  
            IF ( zmask > 0.0 ) vmask_crs(ji,2,jk) = 1.0_wp
               
            zmask = 0.0
            zmask = SUM(umask(ijie,ij:je_2,jk))   
            IF ( zmask > 0.0 ) umask_crs(ji,2,jk) = 1.0_wp
               
            fmask_crs(ji,je_2,jk) = fmask(ijie,2,jk)
         ENDDO
      ENDDO
      !
      DO jk = 1, jpkm1
         DO ji = 2, Nie0_crs  
            ijis = mis_crs(ji)     ;   ijie = mie_crs(ji)       
            DO jj = 3, Nje0_crs
               ijjs = mjs_crs(jj)  ;   ijje = mje_crs(jj)
                          
               zmask = 0.0
               zmask = SUM( tmask(ijis:ijie,ijjs:ijje,jk) ) 
               IF ( zmask > 0.0 ) tmask_crs(ji,jj,jk) = 1.0_wp
               
               zmask = 0.0
               zmask = SUM( vmask(ijis:ijie,ijje     ,jk) )  
               IF ( zmask > 0.0 ) vmask_crs(ji,jj,jk) = 1.0_wp
               
               zmask = 0.0
               zmask = SUM( umask(ijie     ,ijjs:ijje,jk) )  
               IF ( zmask > 0.0 ) umask_crs(ji,jj,jk) = 1.0_wp
               
               fmask_crs(ji,jj,jk) = fmask(ijie,ijje,jk)  
            ENDDO
         ENDDO
      ENDDO

      !
      CALL crs_lbc_lnk( tmask_crs, 'T', 1.0_wp )
      CALL crs_lbc_lnk( vmask_crs, 'V', 1.0_wp )
      CALL crs_lbc_lnk( umask_crs, 'U', 1.0_wp )
      CALL crs_lbc_lnk( fmask_crs, 'F', 1.0_wp )
      !
   END SUBROUTINE crs_dom_msk


   SUBROUTINE crs_dom_coordinates( p_gphi, p_glam, cd_type, p_gphi_crs, p_glam_crs )
      !!----------------------------------------------------------------
      !!               *** SUBROUTINE crs_coordinates ***
      !! ** Purpose :  Determine the coordinates for the coarse grid
      !! 
      !! ** Method  :  From the parent grid subset, search for the central
      !!               point.  For an odd-numbered reduction factor,
      !!               the coordinate will be that of the central T-cell.
      !!               For an even-numbered reduction factor, of a non-square
      !!               coarse grid box, the coordinate will be that of 
      !!               the east or north face or more likely.  For a square
      !!               coarse grid box, the coordinate will be that of
      !!               the central f-corner.
      !!
      !! ** Input   :  p_gphi = parent grid gphi[t|u|v|f]
      !!               p_glam = parent grid glam[t|u|v|f]
      !!               cd_type  = grid type (T,U,V,F)
      !! ** Output  :  p_gphi_crs = coarse grid gphi[t|u|v|f]
      !!               p_glam_crs = coarse grid glam[t|u|v|f]
      !!              
      !! History. 1 Jun.
      !!----------------------------------------------------------------
      !! Arguments
      REAL(wp), DIMENSION(jpi,jpj)        , INTENT(in)  :: p_gphi  ! Parent grid latitude
      REAL(wp), DIMENSION(jpi,jpj)        , INTENT(in)  :: p_glam  ! Parent grid longitude
      CHARACTER(len=1),                     INTENT(in)  :: cd_type   ! grid type (T,U,V,F) 
      REAL(wp), DIMENSION(jpi_crs,jpj_crs), INTENT(out) :: p_gphi_crs  ! Coarse grid latitude
      REAL(wp), DIMENSION(jpi_crs,jpj_crs), INTENT(out) :: p_glam_crs  ! Coarse grid longitude

      !! Local variables
      INTEGER :: ji, jj, jk                   ! dummy loop indices
      INTEGER :: ijis, ijjs

  
      SELECT CASE ( cd_type )
         CASE ( 'T' )
            DO jj =  Njs0_crs, Nje0_crs
               ijjs = mjs_crs(jj) + mybinctr
               DO ji = 2, Nie0_crs
                  ijis = mis_crs(ji) + mxbinctr 
                  p_gphi_crs(ji,jj) = p_gphi(ijis,ijjs)
                  p_glam_crs(ji,jj) = p_glam(ijis,ijjs)
               ENDDO
            ENDDO
         CASE ( 'U' )
            DO jj =  Njs0_crs, Nje0_crs
               ijjs = mjs_crs(jj) + mybinctr                  
               DO ji = 2, Nie0_crs
                  ijis = mis_crs(ji)
                  p_gphi_crs(ji,jj) = p_gphi(ijis,ijjs)
                  p_glam_crs(ji,jj) = p_glam(ijis,ijjs)
               ENDDO
            ENDDO
         CASE ( 'V' )
            DO jj =  Njs0_crs, Nje0_crs
               ijjs = mjs_crs(jj)
               DO ji = 2, Nie0_crs
                  ijis = mis_crs(ji) + mxbinctr 
                  p_gphi_crs(ji,jj) = p_gphi(ijis,ijjs)
                  p_glam_crs(ji,jj) = p_glam(ijis,ijjs)
               ENDDO
            ENDDO
         CASE ( 'F' )
            DO jj =  Njs0_crs, Nje0_crs
               ijjs = mjs_crs(jj)
               DO ji = 2, Nie0_crs
                  ijis = mis_crs(ji)
                  p_gphi_crs(ji,jj) = p_gphi(ijis,ijjs)
                  p_glam_crs(ji,jj) = p_glam(ijis,ijjs)
               ENDDO
            ENDDO
      END SELECT

      ! Retroactively add back the boundary halo cells.
      CALL crs_lbc_lnk( p_gphi_crs, cd_type, 1.0_wp )
      CALL crs_lbc_lnk( p_glam_crs, cd_type, 1.0_wp )
         
      ! Fill up jrow=1 which is zeroed out or not handled by lbc_lnk and lbc_nfd
      SELECT CASE ( cd_type )
         CASE ( 'T', 'V' )
            DO ji = 2, Nie0_crs
               ijis = mis_crs(ji) + mxbinctr 
               p_gphi_crs(ji,1) = p_gphi(ijis,1)
               p_glam_crs(ji,1) = p_glam(ijis,1)
            ENDDO
         CASE ( 'U', 'F' )
            DO ji = 2, Nie0_crs
               ijis = mis_crs(ji) 
               p_gphi_crs(ji,1) = p_gphi(ijis,1)
               p_glam_crs(ji,1) = p_glam(ijis,1)
            ENDDO
      END SELECT
      !
   END SUBROUTINE crs_dom_coordinates

  SUBROUTINE crs_dom_hgr( p_e1, p_e2, cd_type, p_e1_crs, p_e2_crs )
      !!----------------------------------------------------------------
      !!               *** SUBROUTINE crs_dom_hgr ***
      !!
      !! ** Purpose :  Get coarse grid horizontal scale factors and unmasked fraction
      !!
      !! ** Method  :  For grid types T,U,V,Fthe 2D scale factors of 
      !!               the coarse grid are the sum of the east or north faces of the 
      !!               parent grid subset comprising the coarse grid box.     
      !!               - e1,e2 Scale factors
      !!                 Valid arguments: 
      !! ** Inputs  : p_e1, p_e2  = parent grid e1 or e2 (t,u,v,f)
      !!              cd_type     = grid type (T,U,V,F) for scale factors; for velocities (U or V)
      !! ** Outputs : p_e1_crs, p_e2_crs  = parent grid e1 or e2 (t,u,v,f)
      !!
      !! History.     4 Jun.  Write for WGT and scale factors only
      !!----------------------------------------------------------------
      !! 
      !!  Arguments
      REAL(dp), DIMENSION(jpi,jpj)        , INTENT(in)  :: p_e1     ! Parent grid U,V scale factors (e1)
      REAL(dp), DIMENSION(jpi,jpj)        , INTENT(in)  :: p_e2     ! Parent grid U,V scale factors (e2)
Guillaume Samson's avatar
Guillaume Samson committed
      CHARACTER(len=1)                    , INTENT(in)  :: cd_type  ! grid type U,V 

      REAL(wp), DIMENSION(jpi_crs,jpj_crs), INTENT(out) :: p_e1_crs ! Coarse grid box 2D quantity
      REAL(wp), DIMENSION(jpi_crs,jpj_crs), INTENT(out) :: p_e2_crs ! Coarse grid box 2D quantity

      !! Local variables
      INTEGER :: ji, jj, jk     ! dummy loop indices
      INTEGER :: ijie,ijje,ijrs
  
      !!----------------------------------------------------------------  
      ! Initialize      

      DO jk = 1, jpk    
         DO ji = 2, Nie0_crs
            ijie = mie_crs(ji)
            DO jj = Njs0_crs, Nje0_crs
               ijje = mje_crs(jj)   ;   ijrs =  mje_crs(jj) - mjs_crs(jj)
               ! Only for a factro 3 coarsening
               SELECT CASE ( cd_type )
                   CASE ( 'T' )
                      IF( ijrs == 0 .OR. ijrs == 1 ) THEN
                        ! Si à la frontière sud on a pas assez de maille de la grille mère
                         p_e1_crs(ji,jj) = p_e1(ijie-1,ijje) * nn_factx
                         p_e2_crs(ji,jj) = p_e2(ijie-1,ijje) * nn_facty
                      ELSE
                         p_e1_crs(ji,jj) = p_e1(ijie-1,ijje-1) * nn_factx
                         p_e2_crs(ji,jj) = p_e2(ijie-1,ijje-1) * nn_facty
                      ENDIF
                   CASE ( 'U' )
                      IF( ijrs == 0 .OR. ijrs == 1 ) THEN
                         ! Si à la frontière sud on a pas assez de maille de la grille mère
                         p_e1_crs(ji,jj) = p_e1(ijie,ijje) * nn_factx                            
                         p_e2_crs(ji,jj) = p_e2(ijie,ijje) * nn_facty
                      ELSE
                         p_e1_crs(ji,jj) = p_e1(ijie,ijje-1) * nn_factx
                         p_e2_crs(ji,jj) = p_e2(ijie,ijje-1) * nn_facty
                      ENDIF
                   CASE ( 'V' )
                         p_e1_crs(ji,jj) = p_e1(ijie-1,ijje) * nn_factx                            
                         p_e2_crs(ji,jj) = p_e2(ijie-1,ijje) * nn_facty
                   CASE ( 'F' )
                         p_e1_crs(ji,jj) = p_e1(ijie,ijje) * nn_factx                            
                         p_e2_crs(ji,jj) = p_e2(ijie,ijje) * nn_facty
               END SELECT
            ENDDO
         ENDDO
      ENDDO

      CALL crs_lbc_lnk( p_e1_crs, cd_type, 1.0_wp, pfillval=1.0_wp )
      CALL crs_lbc_lnk( p_e2_crs, cd_type, 1.0_wp, pfillval=1.0_wp )

   END SUBROUTINE crs_dom_hgr


   SUBROUTINE crs_dom_facvol( p_mask, cd_type, p_e1, p_e2, p_e3, p_fld1_crs, p_fld2_crs )
      !!----------------------------------------------------------------
      !!               *** SUBROUTINE crsfun_wgt ***
      !! ** Purpose :  Three applications.
      !!               1) SUM. Get coarse grid horizontal scale factors and unmasked fraction
      !!               2) VOL. Get coarse grid box volumes
      !!               3) WGT. Weighting multiplier for volume-weighted and/or
      !!                       area-weighted averages.
      !!                       Weights (i.e. the denominator) calculated here 
      !!                       to avoid IF-tests and division.
      !! ** Method  :  1) SUM.  For grid types T,U,V,F (and W) the 2D scale factors of 
      !!               the coarse grid are the sum of the east or north faces of the 
      !!               parent grid subset comprising the coarse grid box.  
      !!               The fractions of masked:total surface (3D) on the east, 
      !!               north and top faces is, optionally, also output.
      !!               - Top face area sum
      !!                 Valid arguments: cd_type, cd_op='W', p_pmask, p_e1, p_e2 
      !!               - Top face ocean surface fraction
      !!                 Valid arguments: cd_type, cd_op='W', p_pmask, p_e1, p_e2       
      !!               - e1,e2 Scale factors
      !!                 Valid arguments: 
      !!               2) VOL.  For grid types W and T, the coarse grid box
      !!               volumes are output. Also optionally, the fraction of  
      !!               masked:total volume of the parent grid subset is output (i.e. facvol).
      !!               3) WGT. Based on the grid type, the denominator is pre-determined here to  
      !!               perform area- or volume- weighted averages, 
      !!               to avoid IF-tests and divisions.
      !! ** Inputs  : p_e1, p_e2  = parent grid e1 or e2 (t,u,v,f)
      !!              p_pmask     = parent grid mask (T,U,V,F) 
      !!              cd_type     = grid type (T,U,V,F) for scale factors; for velocities (U or V)
      !!              cd_op       = applied operation (SUM, VOL, WGT)
      !!              p_e3      = (Optional) parent grid vertical level thickness (e3u or e3v)
      !! ** Outputs : p_cfield2d_1 = (Optional) 2D field on coarse grid
      !!              p_cfield2d_2 = (Optional) 2D field on coarse grid
      !!              p_cfield3d_1 = (Optional) 3D field on coarse grid
      !!              p_cfield3d_2 = (Optional) 3D field on coarse grid
      !!
      !! History.     4 Jun.  Write for WGT and scale factors only
      !!----------------------------------------------------------------
      CHARACTER(len=1),                         INTENT(in   ) ::   cd_type    ! grid type U,V 
      REAL(wp), DIMENSION(jpi,jpj,jpk)        , INTENT(in   ) ::   p_mask     ! Parent grid U,V mask
      REAL(dp), DIMENSION(jpi,jpj)            , INTENT(in   ) ::   p_e1       ! Parent grid U,V scale factors (e1)
      REAL(dp), DIMENSION(jpi,jpj)            , INTENT(in   ) ::   p_e2       ! Parent grid U,V scale factors (e2)
Guillaume Samson's avatar
Guillaume Samson committed
      REAL(wp), DIMENSION(jpi,jpj,jpk)        , INTENT(in   ) ::   p_e3       ! Parent grid vertical level thickness (e3u, e3v)
      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(  out) ::   p_fld1_crs ! Coarse grid box 3D quantity 
      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(  out) ::   p_fld2_crs ! Coarse grid box 3D quantity 
      !
      INTEGER  ::   ji, jj, jk , ii, ij, je_2
      REAL(wp) ::   zdAm
      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zvol, zmask      
      !!----------------------------------------------------------------  
      !
      !
      p_fld1_crs(:,:,:) = 0._wp
      p_fld2_crs(:,:,:) = 0._wp

      DO jk = 1, jpk
         zvol(:,:,jk) =  p_e1(:,:) * p_e2(:,:) * p_e3(:,:,jk) 
      END DO

      zmask(:,:,:) = 0._wp
      IF( cd_type == 'W' ) THEN
         zmask(:,:,1) = p_mask(:,:,1) 
         DO jk = 2, jpk
            zmask(:,:,jk) = p_mask(:,:,jk-1) 
         ENDDO
      ELSE
         DO jk = 1, jpk
             zmask(:,:,jk) = p_mask(:,:,jk) 
         ENDDO
      ENDIF

      IF( Njs0_crs == 1 .AND. ( ( mje_crs(2) - mjs_crs(2) ) < 2 ) ) THEN     !!cc bande du sud style ORCA2
         IF( mje_crs(2) - mjs_crs(2) == 1 ) THEN
            je_2 = mje_crs(2)
            DO jk = 1, jpk           
               DO ji = nistr, niend, nn_factx
                  ii   = ( ji - mis_crs(2) ) * rfactx_r + 2                 ! cordinate in parent grid
                  p_fld1_crs(ii,2,jk) =  zvol(ji,je_2  ,jk) + zvol(ji+1,je_2  ,jk) + zvol(ji+2,je_2  ,jk)  &
                     &                 + zvol(ji,je_2-1,jk) + zvol(ji+1,je_2-1,jk) + zvol(ji+2,je_2-1,jk) 
                  !
                  zdAm =  zvol(ji  ,je_2,jk) * zmask(ji  ,je_2,jk)  &
                    &   + zvol(ji+1,je_2,jk) * zmask(ji+1,je_2,jk)  &
                    &   + zvol(ji+2,je_2,jk) * zmask(ji+2,je_2,jk) 
                  ! 
                  p_fld2_crs(ii,2,jk) = zdAm / p_fld1_crs(ii,2,jk) 
               ENDDO
            ENDDO
         ENDIF
      ELSE
         je_2 = mjs_crs(2)
         DO jk = 1, jpk           
            DO ji = nistr, niend, nn_factx
               ii   = ( ji - mis_crs(2) ) * rfactx_r + 2    
               p_fld1_crs(ii,2,jk) =  zvol(ji,je_2  ,jk) + zvol(ji+1,je_2  ,jk) + zvol(ji+2,je_2  ,jk)  &
                   &                + zvol(ji,je_2+1,jk) + zvol(ji+1,je_2+1,jk) + zvol(ji+2,je_2+1,jk)  &
                   &                + zvol(ji,je_2+2,jk) + zvol(ji+1,je_2+2,jk) + zvol(ji+2,je_2+2,jk)  
              !
               zdAm = zvol(ji  ,je_2  ,jk) * zmask(ji  ,je_2  ,jk)  &
                 &  + zvol(ji+1,je_2  ,jk) * zmask(ji+1,je_2  ,jk)  &
                 &  + zvol(ji+2,je_2  ,jk) * zmask(ji+2,je_2  ,jk)  &
                 &  + zvol(ji  ,je_2+1,jk) * zmask(ji  ,je_2+1,jk)  &
                 &  + zvol(ji+1,je_2+1,jk) * zmask(ji+1,je_2+1,jk)  &
                 &  + zvol(ji+2,je_2+1,jk) * zmask(ji+2,je_2+1,jk)  &
                 &  + zvol(ji  ,je_2+2,jk) * zmask(ji  ,je_2+2,jk)  &
                 &  + zvol(ji+1,je_2+2,jk) * zmask(ji+1,je_2+2,jk)  &
                 &  + zvol(ji+2,je_2+2,jk) * zmask(ji+2,je_2+2,jk)
                 ! 
                 p_fld2_crs(ii,2,jk) = zdAm / p_fld1_crs(ii,2,jk) 
            ENDDO
         ENDDO
      ENDIF

      DO jk = 1, jpk           
         DO jj  = njstr, njend, nn_facty
            DO ji = nistr, niend, nn_factx
               ii  = ( ji - mis_crs(2) ) * rfactx_r + 2                 ! cordinate in parent grid
               ij  = ( jj - njstr ) * rfacty_r + 3
               !
               p_fld1_crs(ii,ij,jk) =  zvol(ji,jj  ,jk) + zvol(ji+1,jj  ,jk) + zvol(ji+2,jj  ,jk)  &
                   &                 + zvol(ji,jj+1,jk) + zvol(ji+1,jj+1,jk) + zvol(ji+2,jj+1,jk)  &
                   &                 + zvol(ji,jj+2,jk) + zvol(ji+1,jj+2,jk) + zvol(ji+2,jj+2,jk) 
               !
               zdAm =  zvol(ji  ,jj  ,jk) * zmask(ji  ,jj  ,jk)  &
                 &   + zvol(ji+1,jj  ,jk) * zmask(ji+1,jj  ,jk)  &
                 &   + zvol(ji+2,jj  ,jk) * zmask(ji+2,jj  ,jk)  &
                 &   + zvol(ji  ,jj+1,jk) * zmask(ji  ,jj+1,jk)  &
                 &   + zvol(ji+1,jj+1,jk) * zmask(ji+1,jj+1,jk)  &
                 &   + zvol(ji+2,jj+1,jk) * zmask(ji+2,jj+1,jk)  &
                 &   + zvol(ji  ,jj+2,jk) * zmask(ji  ,jj+2,jk)  &
                 &   + zvol(ji+1,jj+2,jk) * zmask(ji+1,jj+2,jk)  &
                 &   + zvol(ji+2,jj+2,jk) * zmask(ji+2,jj+2,jk)
                 ! 
                p_fld2_crs(ii,ij,jk) = zdAm / p_fld1_crs(ii,ij,jk) 
            ENDDO
         ENDDO
      ENDDO
      !                                             !  Retroactively add back the boundary halo cells.
      CALL crs_lbc_lnk( p_fld1_crs, cd_type, 1.0_wp ) 
      CALL crs_lbc_lnk( p_fld2_crs, cd_type, 1.0_wp ) 
      !
      !
   END SUBROUTINE crs_dom_facvol


   SUBROUTINE crs_dom_ope_3d( p_fld, cd_op, cd_type, p_mask, p_fld_crs, p_e12, p_e3, p_surf_crs, p_mask_crs, psgn )
      !!----------------------------------------------------------------
      !!               *** SUBROUTINE crsfun_UV ***
      !! ** Purpose :  Average, area-weighted, of U or V on the east and north faces 
      !!
      !! ** Method  :  The U and V velocities (3D) are determined as the area-weighted averages
      !!               on the east and north faces, respectively,
      !!               of the parent grid subset comprising the coarse grid box. 
      !!               In the case of the V and F grid, the last jrow minus 1 is spurious.
      !! ** Inputs  : p_e1_e2     = parent grid e1 or e2 (t,u,v,f)
      !!              cd_type     = grid type (T,U,V,F) for scale factors; for velocities (U or V)
      !!              psgn        = sign change over north fold (See lbclnk.F90)
      !!              p_pmask     = parent grid mask (T,U,V,F) for scale factors; 
      !!                                       for velocities (U or V)
      !!              p_e3      = parent grid vertical level thickness (e3u or e3v)
      !!              p_pfield    = U or V on the parent grid
      !!              p_surf_crs  = (Optional) Coarse grid weight for averaging
      !! ** Outputs : p_cfield3d = 3D field on coarse grid
      !!
      !! History.  29 May.  completed draft.
      !!            4 Jun.  Revision for WGT
      !!            5 Jun.  Streamline for area-weighted average only ; separate scale factor and weights.
      !!----------------------------------------------------------------
      REAL(wp), DIMENSION(jpi,jpj,jpk),         INTENT(in)           :: p_fld      ! T, U, V or W on parent grid
      CHARACTER(len=3),                         INTENT(in)           :: cd_op      ! Operation SUM, MAX or MIN
      CHARACTER(len=1),                         INTENT(in)           :: cd_type    ! grid type U,V 
      REAL(wp), DIMENSION(jpi,jpj,jpk),         INTENT(in)           :: p_mask     ! Parent grid T,U,V mask
      REAL(dp), DIMENSION(jpi,jpj),             INTENT(in), OPTIONAL :: p_e12      ! Parent grid T,U,V scale factors (e1 or e2)
Guillaume Samson's avatar
Guillaume Samson committed
477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000
      REAL(wp), DIMENSION(jpi,jpj,jpk),         INTENT(in), OPTIONAL :: p_e3       ! Parent grid vertical level thickness (e3u, e3v)
      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(in), OPTIONAL :: p_surf_crs ! Coarse grid area-weighting denominator    
      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(in), OPTIONAL :: p_mask_crs ! Coarse grid T,U,V maska
      REAL(wp),                                 INTENT(in)           :: psgn       ! sign 
      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(  out)        :: p_fld_crs  ! Coarse grid box 3D quantity 
      !
      INTEGER  ::   ji, jj, jk 
      INTEGER  ::   ii, ij, ijie, ijje, je_2
      REAL(wp) ::   zflcrs, zsfcrs   
      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zsurf, zsurfmsk, zmask  
      !!----------------------------------------------------------------  
      !
      p_fld_crs(:,:,:) = 0._wp
      !
      SELECT CASE ( cd_op )
      !
      CASE ( 'VOL' )
         !
         ALLOCATE( zsurf(jpi,jpj,jpk), zsurfmsk(jpi,jpj,jpk) )
         !   
         SELECT CASE ( cd_type )
         !   
         CASE( 'T', 'W' )
               IF( cd_type == 'T' ) THEN
                     DO jk = 1, jpk
                        zsurf   (:,:,jk) =  p_e12(:,:) * p_e3(:,:,jk) *  p_mask(:,:,jk) 
                        zsurfmsk(:,:,jk) =  zsurf(:,:,jk) 
                    ENDDO
                  ELSE
                     zsurf   (:,:,1) =  p_e12(:,:) * p_e3(:,:,1)
                     zsurfmsk(:,:,1) =  zsurf(:,:,1) *  p_mask(:,:,1) 
                     DO jk = 2, jpk
                        zsurf   (:,:,jk) =  p_e12(:,:) * p_e3(:,:,jk)
                        zsurfmsk(:,:,jk) =  zsurf(:,:,jk) * p_mask(:,:,jk-1) 
                     ENDDO
                  ENDIF
         
                  IF( Njs0_crs == 1 .AND. ( ( mje_crs(2) - mjs_crs(2) ) < 2 ) ) THEN     !!cc bande du sud style ORCA2
                     IF( mje_crs(2) - mjs_crs(2) == 1 ) THEN
                        je_2 = mje_crs(2)
                        DO jk = 1, jpk           
                           DO ji = nistr, niend, nn_factx
                              ii   = ( ji - mis_crs(2) ) * rfactx_r + 2       
                              zflcrs =  p_fld(ji  ,je_2,jk) * zsurfmsk(ji  ,je_2,jk)   &
                                &     + p_fld(ji+1,je_2,jk) * zsurfmsk(ji+1,je_2,jk)   &
                                &     + p_fld(ji+2,je_2,jk) * zsurfmsk(ji+2,je_2,jk) 
 
                              zsfcrs =  zsurf(ji,je_2,jk) + zsurf(ji+1,je_2,jk) + zsurf(ji+2,je_2,jk) 
                              !
                              p_fld_crs(ii,2,jk) = zflcrs
                              IF( zsfcrs /= 0.0 )  p_fld_crs(ii,2,jk) = zflcrs / zsfcrs
                           ENDDO
                        ENDDO
                     ENDIF
                  ELSE
                     je_2 = mjs_crs(2)
                     DO jk = 1, jpk           
                        DO ji = nistr, niend, nn_factx
                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2   
                           zflcrs =  p_fld(ji  ,je_2  ,jk) * zsurfmsk(ji  ,je_2  ,jk) &
                             &     + p_fld(ji+1,je_2  ,jk) * zsurfmsk(ji+1,je_2  ,jk) &
                             &     + p_fld(ji+2,je_2  ,jk) * zsurfmsk(ji+2,je_2  ,jk) &
                             &     + p_fld(ji  ,je_2+1,jk) * zsurfmsk(ji  ,je_2+1,jk) &
                             &     + p_fld(ji+1,je_2+1,jk) * zsurfmsk(ji+1,je_2+1,jk) &
                             &     + p_fld(ji+2,je_2+1,jk) * zsurfmsk(ji+2,je_2+1,jk) &
                             &     + p_fld(ji  ,je_2+2,jk) * zsurfmsk(ji  ,je_2+2,jk) &
                             &     + p_fld(ji+1,je_2+2,jk) * zsurfmsk(ji+1,je_2+2,jk) &
                             &     + p_fld(ji+2,je_2+2,jk) * zsurfmsk(ji+2,je_2+2,jk) 

                           zsfcrs =  zsurf(ji,je_2  ,jk) + zsurf(ji+1,je_2  ,jk) + zsurf(ji+2,je_2  ,jk) &
                             &     + zsurf(ji,je_2+1,jk) + zsurf(ji+1,je_2+1,jk) + zsurf(ji+2,je_2+1,jk) &
                             &     + zsurf(ji,je_2+2,jk) + zsurf(ji+1,je_2+2,jk) + zsurf(ji+2,je_2+2,jk) 
                            !
                            p_fld_crs(ii,2,jk) = zflcrs
                            IF( zsfcrs /= 0.0 )  p_fld_crs(ii,2,jk) = zflcrs / zsfcrs
                        ENDDO
                     ENDDO
                  ENDIF
                  !
                  DO jk = 1, jpk           
                     DO jj  = njstr, njend, nn_facty
                        DO ji = nistr, niend, nn_factx
                           ii = ( ji - mis_crs(2) ) * rfactx_r + 2                 ! cordinate in parent grid
                           ij = ( jj - njstr ) * rfacty_r + 3
                           zflcrs =  p_fld(ji  ,jj  ,jk) * zsurfmsk(ji  ,jj  ,jk) &
                             &     + p_fld(ji+1,jj  ,jk) * zsurfmsk(ji+1,jj  ,jk) &
                             &     + p_fld(ji+2,jj  ,jk) * zsurfmsk(ji+2,jj  ,jk) &
                             &     + p_fld(ji  ,jj+1,jk) * zsurfmsk(ji  ,jj+1,jk) &
                             &     + p_fld(ji+1,jj+1,jk) * zsurfmsk(ji+1,jj+1,jk) &
                             &     + p_fld(ji+2,jj+1,jk) * zsurfmsk(ji+2,jj+1,jk) &
                             &     + p_fld(ji  ,jj+2,jk) * zsurfmsk(ji  ,jj+2,jk) &
                             &     + p_fld(ji+1,jj+2,jk) * zsurfmsk(ji+1,jj+2,jk) &
                             &     + p_fld(ji+2,jj+2,jk) * zsurfmsk(ji+2,jj+2,jk) 

                           zsfcrs =  zsurf(ji,jj  ,jk) + zsurf(ji+1,jj  ,jk) + zsurf(ji+2,jj  ,jk) &
                             &     + zsurf(ji,jj+1,jk) + zsurf(ji+1,jj+1,jk) + zsurf(ji+2,jj+1,jk) &
                             &     + zsurf(ji,jj+2,jk) + zsurf(ji+1,jj+2,jk) + zsurf(ji+2,jj+2,jk) 
                            !
                           p_fld_crs(ii,ij,jk) = zflcrs
                           IF( zsfcrs /= 0.0 )  p_fld_crs(ii,ij,jk) = zflcrs / zsfcrs
                        ENDDO      
                     ENDDO
                  ENDDO  
               CASE DEFAULT
                  CALL ctl_stop( 'STOP', 'error from crs_dom_ope_3d, you should not be there...' )
               END SELECT

              DEALLOCATE( zsurf, zsurfmsk )

         CASE ( 'SUM' )
         
            ALLOCATE( zsurfmsk(jpi,jpj,jpk) )

            SELECT CASE ( cd_type )
              CASE( 'W' )
                  IF( PRESENT( p_e3 ) ) THEN
                    zsurfmsk(:,:,1) =  p_e12(:,:) * p_e3(:,:,1) * p_mask(:,:,1) 
                    DO jk = 2, jpk
                      zsurfmsk(:,:,jk) =  p_e12(:,:) * p_e3(:,:,jk) * p_mask(:,:,jk-1) 
                    ENDDO
                 ELSE
                    zsurfmsk(:,:,1) =  p_e12(:,:) * p_mask(:,:,1) 
                    DO jk = 2, jpk
                      zsurfmsk(:,:,jk) =  p_e12(:,:) * p_mask(:,:,jk-1) 
                    ENDDO
                 ENDIF
              CASE DEFAULT
                 IF( PRESENT( p_e3 ) ) THEN
                    DO jk = 1, jpk
                      zsurfmsk(:,:,jk) =  p_e12(:,:) * p_e3(:,:,jk) * p_mask(:,:,jk) 
                    ENDDO
                 ELSE
                    DO jk = 1, jpk
                      zsurfmsk(:,:,jk) =  p_e12(:,:) * p_mask(:,:,jk) 
                    ENDDO
                 ENDIF
              END SELECT

            SELECT CASE ( cd_type )
            
               CASE( 'T', 'W' )
         
                  IF( Njs0_crs == 1 .AND. ( ( mje_crs(2) - mjs_crs(2) ) < 2 ) ) THEN     !!cc bande du sud style ORCA2
                     IF( mje_crs(2) - mjs_crs(2) == 1 ) THEN
                        je_2 = mje_crs(2)
                        DO jk = 1, jpk           
                           DO ji = nistr, niend, nn_factx
                              ii   = ( ji - mis_crs(2) ) * rfactx_r + 2          
                              zflcrs  =  p_fld(ji  ,je_2,jk) * zsurfmsk(ji  ,je_2,jk) &
                                &      + p_fld(ji+1,je_2,jk) * zsurfmsk(ji+1,je_2,jk) &
                                &      + p_fld(ji+2,je_2,jk) * zsurfmsk(ji+2,je_2,jk) 
                               !
                              p_fld_crs(ii,2,jk) = zflcrs
                           ENDDO
                        ENDDO
                      ENDIF
                  ELSE
                     je_2 = mjs_crs(2)
                     DO jk = 1, jpk           
                        DO ji = nistr, niend, nn_factx
                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2   
                           zflcrs  =  p_fld(ji  ,je_2  ,jk) * zsurfmsk(ji  ,je_2  ,jk)  &
                             &      + p_fld(ji+1,je_2  ,jk) * zsurfmsk(ji+1,je_2  ,jk)  &
                             &      + p_fld(ji+2,je_2  ,jk) * zsurfmsk(ji+2,je_2  ,jk)  &
                             &      + p_fld(ji  ,je_2+1,jk) * zsurfmsk(ji  ,je_2+1,jk)  &
                             &      + p_fld(ji+1,je_2+1,jk) * zsurfmsk(ji+1,je_2+1,jk)  &
                             &      + p_fld(ji+2,je_2+1,jk) * zsurfmsk(ji+2,je_2+1,jk)  &
                             &      + p_fld(ji  ,je_2+2,jk) * zsurfmsk(ji  ,je_2+2,jk)  &
                             &      + p_fld(ji+1,je_2+2,jk) * zsurfmsk(ji+1,je_2+2,jk)  &
                             &      + p_fld(ji+2,je_2+2,jk) * zsurfmsk(ji+2,je_2+2,jk)  
                            !
                            p_fld_crs(ii,2,jk) = zflcrs
                        ENDDO
                     ENDDO
                  ENDIF
                  !
                  DO jk = 1, jpk           
                     DO jj  = njstr, njend, nn_facty
                        DO ji = nistr, niend, nn_factx
                           ii  = ( ji - mis_crs(2) ) * rfactx_r + 2                 ! cordinate in parent grid
                           ij  = ( jj - njstr ) * rfacty_r + 3
                           zflcrs  =  p_fld(ji  ,jj  ,jk) * zsurfmsk(ji  ,jj  ,jk)  &
                             &      + p_fld(ji+1,jj  ,jk) * zsurfmsk(ji+1,jj  ,jk)  &
                             &      + p_fld(ji+2,jj  ,jk) * zsurfmsk(ji+2,jj  ,jk)  &
                             &      + p_fld(ji  ,jj+1,jk) * zsurfmsk(ji  ,jj+1,jk)  &
                             &      + p_fld(ji+1,jj+1,jk) * zsurfmsk(ji+1,jj+1,jk)  &
                             &      + p_fld(ji+2,jj+1,jk) * zsurfmsk(ji+2,jj+1,jk)  &
                             &      + p_fld(ji  ,jj+2,jk) * zsurfmsk(ji  ,jj+2,jk)  &
                             &      + p_fld(ji+1,jj+2,jk) * zsurfmsk(ji+1,jj+2,jk)  &
                             &      + p_fld(ji+2,jj+2,jk) * zsurfmsk(ji+2,jj+2,jk)  
                            !
                            p_fld_crs(ii,ij,jk) = zflcrs
                            ! 
                        ENDDO      
                     ENDDO
                  ENDDO   
            
               CASE( 'V' )

                  IF( Njs0_crs == 1 .AND. ( ( mje_crs(2) - mjs_crs(2) ) < 2 ) ) THEN     !!cc bande du sud style ORCA2
                     IF( mje_crs(2) - mjs_crs(2) == 1 ) THEN
                        ijje = mje_crs(2)
                     ENDIF
                  ELSE
                     ijje = mjs_crs(2)
                  ENDIF
                  !
                  DO jk = 1, jpk           
                     DO ji = nistr, niend, nn_factx
                        ii   = ( ji - mis_crs(2) ) * rfactx_r + 2    
                        zflcrs  =  p_fld(ji  ,ijje,jk) * zsurfmsk(ji  ,ijje,jk) &
                          &      + p_fld(ji+1,ijje,jk) * zsurfmsk(ji+1,ijje,jk) &
                          &      + p_fld(ji+2,ijje,jk) * zsurfmsk(ji+2,ijje,jk) 
                          !
                        p_fld_crs(ii,2,jk) = zflcrs
                     ENDDO
                  ENDDO
                  !
                  DO jk = 1, jpk           
                     DO jj  = njstr, njend, nn_facty
                        DO ji = nistr, niend, nn_factx
                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2                 ! cordinate in parent grid
                           ij   = ( jj - njstr ) * rfacty_r + 3
                           ijje = mje_crs(ij)
                           zflcrs  =  p_fld(ji  ,ijje,jk) * zsurfmsk(ji  ,ijje,jk) &
                             &      + p_fld(ji+1,ijje,jk) * zsurfmsk(ji+1,ijje,jk) &
                             &      + p_fld(ji+2,ijje,jk) * zsurfmsk(ji+2,ijje,jk) 
                             !
                           p_fld_crs(ii,ij,jk) = zflcrs
                           ! 
                        ENDDO      
                     ENDDO
                  ENDDO   
            
               CASE( 'U' )

                  IF( Njs0_crs == 1 .AND. ( ( mje_crs(2) - mjs_crs(2) ) < 2 ) ) THEN     !!cc bande du sud style ORCA2
                     IF( mje_crs(2) - mjs_crs(2) == 1 ) THEN
                        je_2 = mje_crs(2)
                        DO jk = 1, jpk           
                           DO ji = nistr, niend, nn_factx
                              ii   = ( ji - mis_crs(2) ) * rfactx_r + 2  
                              ijie = mie_crs(ii)
                              zflcrs  =  p_fld(ijie,je_2,jk) * zsurfmsk(ijie,je_2,jk)  
                              p_fld_crs(ii,2,jk) = zflcrs
                           ENDDO
                        ENDDO
                      ENDIF
                  ELSE
                     je_2 = mjs_crs(2)
                     DO jk = 1, jpk           
                        DO ji = nistr, niend, nn_factx
                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2   
                           ijie = mie_crs(ii)
                           zflcrs =  p_fld(ijie,je_2  ,jk) * zsurfmsk(ijie,je_2  ,jk)  &
                             &     + p_fld(ijie,je_2+1,jk) * zsurfmsk(ijie,je_2+1,jk)  &
                             &     + p_fld(ijie,je_2+2,jk) * zsurfmsk(ijie,je_2+2,jk) 

                           p_fld_crs(ii,2,jk) = zflcrs
                        ENDDO
                     ENDDO
                  ENDIF
                  !
                  DO jk = 1, jpk           
                     DO jj  = njstr, njend, nn_facty
                        DO ji = nistr, niend, nn_factx
                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
                           ij   = ( jj - njstr ) * rfacty_r + 3
                           ijie = mie_crs(ii)
                           zflcrs =  p_fld(ijie,jj  ,jk) * zsurfmsk(ijie,jj  ,jk)  &
                              &    + p_fld(ijie,jj+1,jk) * zsurfmsk(ijie,jj+1,jk)  &
                              &    + p_fld(ijie,jj+2,jk) * zsurfmsk(ijie,jj+2,jk) 
                             !
                           p_fld_crs(ii,ij,jk) = zflcrs
                           ! 
                        ENDDO      
                     ENDDO
                  ENDDO   

              END SELECT

              IF( PRESENT( p_surf_crs ) ) THEN
                 WHERE ( p_surf_crs /= 0.0 ) p_fld_crs(:,:,:) = p_fld_crs(:,:,:) / p_surf_crs(:,:,:)
              ENDIF

              DEALLOCATE( zsurfmsk )

         CASE ( 'MAX' )    !  search the max of unmasked grid cells

            ALLOCATE( zmask(jpi,jpj,jpk) )

            SELECT CASE ( cd_type )
              CASE( 'W' )
                  zmask(:,:,1) = p_mask(:,:,1) 
                  DO jk = 2, jpk
                     zmask(:,:,jk) = p_mask(:,:,jk-1) 
                  ENDDO
              CASE ( 'T' )
                  DO jk = 1, jpk
                     zmask(:,:,jk) = p_mask(:,:,jk) 
                  ENDDO
            END SELECT

            SELECT CASE ( cd_type )
            
               CASE( 'T', 'W' )
         
                  IF( Njs0_crs == 1 .AND. ( ( mje_crs(2) - mjs_crs(2) ) < 2 ) ) THEN     !!cc bande du sud style ORCA2
                     IF( mje_crs(2) - mjs_crs(2) == 1 ) THEN
                        je_2 = mje_crs(2)
                        DO jk = 1, jpk           
                           DO ji = nistr, niend, nn_factx
                              ii   = ( ji - mis_crs(2) ) * rfactx_r + 2        
                              zflcrs =  &
                                & MAX( p_fld(ji  ,je_2,jk) * zmask(ji  ,je_2,jk) - ( 1.- zmask(ji  ,je_2,jk) ) * r_inf ,  &
                                &      p_fld(ji+1,je_2,jk) * zmask(ji+1,je_2,jk) - ( 1.- zmask(ji+1,je_2,jk) ) * r_inf ,  &
                                &      p_fld(ji+2,je_2,jk) * zmask(ji+2,je_2,jk) - ( 1.- zmask(ji+2,je_2,jk) ) * r_inf  )
                              !
                              p_fld_crs(ii,2,jk) = zflcrs
                           ENDDO
                        ENDDO
                      ENDIF
                  ELSE
                     je_2 = mjs_crs(2)
                     DO jk = 1, jpk           
                        DO ji = nistr, niend, nn_factx
                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2   
                           zflcrs =  &
                             & MAX( p_fld(ji  ,je_2  ,jk) * zmask(ji  ,je_2  ,jk) - ( 1.- zmask(ji  ,je_2  ,jk) ) * r_inf ,  &
                             &      p_fld(ji+1,je_2  ,jk) * zmask(ji+1,je_2  ,jk) - ( 1.- zmask(ji+1,je_2  ,jk) ) * r_inf ,  &
                             &      p_fld(ji+2,je_2  ,jk) * zmask(ji+2,je_2  ,jk) - ( 1.- zmask(ji+2,je_2  ,jk) ) * r_inf ,  &
                             &      p_fld(ji  ,je_2+1,jk) * zmask(ji  ,je_2+1,jk) - ( 1.- zmask(ji  ,je_2+1,jk) ) * r_inf ,  &
                             &      p_fld(ji+1,je_2+1,jk) * zmask(ji+1,je_2+1,jk) - ( 1.- zmask(ji+1,je_2+1,jk) ) * r_inf ,  &
                             &      p_fld(ji+2,je_2+1,jk) * zmask(ji+2,je_2+1,jk) - ( 1.- zmask(ji+2,je_2+1,jk) ) * r_inf ,  &
                             &      p_fld(ji  ,je_2+2,jk) * zmask(ji  ,je_2+2,jk) - ( 1.- zmask(ji  ,je_2+2,jk) ) * r_inf ,  &
                             &      p_fld(ji+1,je_2+2,jk) * zmask(ji+1,je_2+2,jk) - ( 1.- zmask(ji+1,je_2+2,jk) ) * r_inf ,  &
                             &      p_fld(ji+2,je_2+2,jk) * zmask(ji+2,je_2+2,jk) - ( 1.- zmask(ji+2,je_2+2,jk) ) * r_inf   )
                           !
                           p_fld_crs(ii,2,jk) = zflcrs
                        ENDDO
                     ENDDO
                  ENDIF
                  !
                  DO jk = 1, jpk           
                     DO jj  = njstr, njend, nn_facty
                        DO ji = nistr, niend, nn_factx
                           ii  = ( ji - mis_crs(2) ) * rfactx_r + 2                 ! cordinate in parent grid
                           ij  = ( jj - njstr ) * rfacty_r + 3
                           zflcrs =  &
                             & MAX( p_fld(ji  ,jj  ,jk) * zmask(ji  ,jj  ,jk) - ( 1.- zmask(ji  ,jj  ,jk) ) * r_inf ,  &
                             &      p_fld(ji+1,jj  ,jk) * zmask(ji+1,jj  ,jk) - ( 1.- zmask(ji+1,jj  ,jk) ) * r_inf ,  &
                             &      p_fld(ji+2,jj  ,jk) * zmask(ji+2,jj  ,jk) - ( 1.- zmask(ji+2,jj  ,jk) ) * r_inf ,  &
                             &      p_fld(ji  ,jj+1,jk) * zmask(ji  ,jj+1,jk) - ( 1.- zmask(ji  ,jj+1,jk) ) * r_inf ,  &
                             &      p_fld(ji+1,jj+1,jk) * zmask(ji+1,jj+1,jk) - ( 1.- zmask(ji+1,jj+1,jk) ) * r_inf ,  &
                             &      p_fld(ji+2,jj+1,jk) * zmask(ji+2,jj+1,jk) - ( 1.- zmask(ji+2,jj+1,jk) ) * r_inf ,  &
                             &      p_fld(ji  ,jj+2,jk) * zmask(ji  ,jj+2,jk) - ( 1.- zmask(ji  ,jj+2,jk) ) * r_inf ,  &
                             &      p_fld(ji+1,jj+2,jk) * zmask(ji+1,jj+2,jk) - ( 1.- zmask(ji+1,jj+2,jk) ) * r_inf ,  &
                             &      p_fld(ji+2,jj+2,jk) * zmask(ji+2,jj+2,jk) - ( 1.- zmask(ji+2,jj+2,jk) ) * r_inf   )
                           !
                           p_fld_crs(ii,ij,jk) = zflcrs
                           !
                        ENDDO      
                     ENDDO
                  ENDDO   
            
               CASE( 'V' )

                 IF( Njs0_crs == 1 .AND. ( ( mje_crs(2) - mjs_crs(2) ) < 2 ) ) THEN     !!cc bande du sud style ORCA2
                     IF( mje_crs(2) - mjs_crs(2) == 1 ) THEN
                        ijje = mje_crs(2)
                      ENDIF
                  ELSE
                     ijje = mjs_crs(2)
                  ENDIF

                  DO jk = 1, jpk
                     DO ji = nistr, niend, nn_factx
                        ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
                        zflcrs = &
                          & MAX( p_fld(ji  ,ijje,jk) * p_mask(ji  ,ijje,jk) - ( 1.- p_mask(ji,ijje,jk) ) * r_inf ,  &
                          &      p_fld(ji+1,ijje,jk) * p_mask(ji+1,ijje,jk) - ( 1.- p_mask(ji,ijje,jk) ) * r_inf ,  &
                          &      p_fld(ji+2,ijje,jk) * p_mask(ji+2,ijje,jk) - ( 1.- p_mask(ji,ijje,jk) ) * r_inf )
                          !
                        p_fld_crs(ii,2,jk) = zflcrs
                     ENDDO
                  ENDDO
                  !
                  DO jk = 1, jpk           
                     DO jj  = njstr, njend, nn_facty
                        DO ji = nistr, niend, nn_factx
                           ii  = ( ji - mis_crs(2) ) * rfactx_r + 2                 ! cordinate in parent grid
                           ij  = ( jj - njstr ) * rfacty_r + 3
                           ijje = mje_crs(ij)
                           !                  
                           zflcrs = &
                             & MAX( p_fld(ji  ,ijje,jk) * p_mask(ji  ,ijje,jk) - ( 1.- p_mask(ji,ijje,jk) ) * r_inf ,  &
                             &      p_fld(ji+1,ijje,jk) * p_mask(ji+1,ijje,jk) - ( 1.- p_mask(ji,ijje,jk) ) * r_inf ,  &
                             &      p_fld(ji+2,ijje,jk) * p_mask(ji+2,ijje,jk) - ( 1.- p_mask(ji,ijje,jk) ) * r_inf )
                           !
                           p_fld_crs(ii,ij,jk) = zflcrs
                           !
                        ENDDO      
                     ENDDO
                  ENDDO   

            
               CASE( 'U' )

                 IF( Njs0_crs == 1 .AND. ( ( mje_crs(2) - mjs_crs(2) ) < 2 ) ) THEN     !!cc bande du sud style ORCA2
                     IF( mje_crs(2) - mjs_crs(2) == 1 ) THEN
                        je_2 = mje_crs(2)
                        DO jk = 1, jpk           
                           DO ji = nistr, niend, nn_factx
                              ii   = ( ji - mis_crs(2) ) * rfactx_r + 2  
                              ijie = mie_crs(ii)
                              zflcrs = p_fld(ijie,je_2,jk) * p_mask(ijie,je_2,jk) - ( 1.- p_mask(ijie,je_2,jk) ) * r_inf 
                              !
                              p_fld_crs(ii,2,jk) = zflcrs
                            ENDDO
                        ENDDO
                      ENDIF
                  ELSE
                     je_2 = mjs_crs(2)
                     DO jk = 1, jpk           
                        DO ji = nistr, niend, nn_factx
                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2   
                           ijie = mie_crs(ii)
                           zflcrs = &
                             & MAX( p_fld(ijie,je_2  ,jk) * p_mask(ijie,je_2  ,jk) - ( 1.- p_mask(ijie,je_2,jk) ) * r_inf ,  &
                             &      p_fld(ijie,je_2+1,jk) * p_mask(ijie,je_2+1,jk) - ( 1.- p_mask(ijie,je_2,jk) ) * r_inf ,  &
                             &      p_fld(ijie,je_2+2,jk) * p_mask(ijie,je_2+2,jk) - ( 1.- p_mask(ijie,je_2,jk) ) * r_inf  )
                            !
                           p_fld_crs(ii,2,jk) = zflcrs
                        ENDDO
                     ENDDO
                  ENDIF
                  !
                  DO jk = 1, jpk           
                     DO jj  = njstr, njend, nn_facty
                        DO ji = nistr, niend, nn_factx
                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2     
                           ij   = ( jj - njstr ) * rfacty_r + 3
                           ijie = mie_crs(ii)
                           zflcrs =  &
                             & MAX( p_fld(ijie,jj  ,jk) * p_mask(ijie,jj  ,jk) - ( 1.- p_mask(ijie,jj,jk) ) * r_inf ,  &
                             &      p_fld(ijie,jj+1,jk) * p_mask(ijie,jj+1,jk) - ( 1.- p_mask(ijie,jj,jk) ) * r_inf ,  &
                             &      p_fld(ijie,jj+2,jk) * p_mask(ijie,jj+2,jk) - ( 1.- p_mask(ijie,jj,jk) ) * r_inf  )
                           ! 
                           p_fld_crs(ii,ij,jk) = zflcrs
                           ! 
                        ENDDO      
                     ENDDO
                  ENDDO   

              END SELECT

              DEALLOCATE( zmask )

         CASE ( 'MIN' )      !   Search the min of unmasked grid cells

            ALLOCATE( zmask(jpi,jpj,jpk) )

            SELECT CASE ( cd_type )
              CASE( 'W' )
                  zmask(:,:,1) = p_mask(:,:,1) 
                  DO jk = 2, jpk
                     zmask(:,:,jk) = p_mask(:,:,jk-1) 
                  ENDDO
              CASE ( 'T' )
                  DO jk = 1, jpk
                     zmask(:,:,jk) = p_mask(:,:,jk) 
                  ENDDO
            END SELECT

            SELECT CASE ( cd_type )

               CASE( 'T', 'W' )
         
                  IF( Njs0_crs == 1 .AND. ( ( mje_crs(2) - mjs_crs(2) ) < 2 ) ) THEN     !!cc bande du sud style ORCA2
                     IF( mje_crs(2) - mjs_crs(2) == 1 ) THEN
                        je_2 = mje_crs(2)
                        DO jk = 1, jpk           
                           DO ji = nistr, niend, nn_factx
                              ii   = ( ji - mis_crs(2) ) * rfactx_r + 2        
                              zflcrs =  &
                                & MIN( p_fld(ji  ,je_2,jk) * zmask(ji  ,je_2,jk) + ( 1.- zmask(ji  ,je_2,jk) ) * r_inf ,  &
                                &      p_fld(ji+1,je_2,jk) * zmask(ji+1,je_2,jk) + ( 1.- zmask(ji+1,je_2,jk) ) * r_inf ,  &
                                &      p_fld(ji+2,je_2,jk) * zmask(ji+2,je_2,jk) + ( 1.- zmask(ji+2,je_2,jk) ) * r_inf  )
                              !
                              p_fld_crs(ii,2,jk) = zflcrs
                           ENDDO
                        ENDDO
                      ENDIF
                  ELSE
                     je_2 = mjs_crs(2)
                     DO jk = 1, jpk           
                        DO ji = nistr, niend, nn_factx
                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2   
                           zflcrs =  &
                             & MIN( p_fld(ji  ,je_2  ,jk) * zmask(ji  ,je_2  ,jk) + ( 1.- zmask(ji  ,je_2  ,jk) ) * r_inf ,  &
                             &      p_fld(ji+1,je_2  ,jk) * zmask(ji+1,je_2  ,jk) + ( 1.- zmask(ji+1,je_2  ,jk) ) * r_inf ,  &
                             &      p_fld(ji+2,je_2  ,jk) * zmask(ji+2,je_2  ,jk) + ( 1.- zmask(ji+2,je_2  ,jk) ) * r_inf ,  &
                             &      p_fld(ji  ,je_2+1,jk) * zmask(ji  ,je_2+1,jk) + ( 1.- zmask(ji  ,je_2+1,jk) ) * r_inf ,  &
                             &      p_fld(ji+1,je_2+1,jk) * zmask(ji+1,je_2+1,jk) + ( 1.- zmask(ji+1,je_2+1,jk) ) * r_inf ,  &
                             &      p_fld(ji+2,je_2+1,jk) * zmask(ji+2,je_2+1,jk) + ( 1.- zmask(ji+2,je_2+1,jk) ) * r_inf ,  &
                             &      p_fld(ji  ,je_2+2,jk) * zmask(ji  ,je_2+2,jk) + ( 1.- zmask(ji  ,je_2+2,jk) ) * r_inf ,  &
                             &      p_fld(ji+1,je_2+2,jk) * zmask(ji+1,je_2+2,jk) + ( 1.- zmask(ji+1,je_2+2,jk) ) * r_inf ,  &
                             &      p_fld(ji+2,je_2+2,jk) * zmask(ji+2,je_2+2,jk) + ( 1.- zmask(ji+2,je_2+2,jk) ) * r_inf   )
                           !
                           p_fld_crs(ii,2,jk) = zflcrs
                        ENDDO
                     ENDDO
                  ENDIF
                  !
                  DO jk = 1, jpk           
                     DO jj  = njstr, njend, nn_facty
                        DO ji = nistr, niend, nn_factx
                           ii  = ( ji - mis_crs(2) ) * rfactx_r + 2                 ! cordinate in parent grid
                           ij  = ( jj - njstr ) * rfacty_r + 3
                           zflcrs =  &
                             & MIN( p_fld(ji  ,jj  ,jk) * zmask(ji  ,jj  ,jk) + ( 1.- zmask(ji  ,jj  ,jk) ) * r_inf ,  &
                             &      p_fld(ji+1,jj  ,jk) * zmask(ji+1,jj  ,jk) + ( 1.- zmask(ji+1,jj  ,jk) ) * r_inf ,  &
                             &      p_fld(ji+2,jj  ,jk) * zmask(ji+2,jj  ,jk) + ( 1.- zmask(ji+2,jj  ,jk) ) * r_inf ,  &
                             &      p_fld(ji  ,jj+1,jk) * zmask(ji  ,jj+1,jk) + ( 1.- zmask(ji  ,jj+1,jk) ) * r_inf ,  &