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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 11573 for NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/BDY/bdylib.F90 – NEMO

Ignore:
Timestamp:
2019-09-19T11:18:03+02:00 (5 years ago)
Author:
jchanut
Message:

#2222, merged with trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/BDY/bdylib.F90

    r10529 r11573  
    1515   USE bdy_oce        ! ocean open boundary conditions 
    1616   USE phycst         ! physical constants 
     17   USE bdyini 
    1718   ! 
    1819   USE in_out_manager ! 
     
    7576      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta  ! tracer trend 
    7677      !! 
    77       REAL(wp) ::   zwgt           ! boundary weight 
    7878      INTEGER  ::   ib, ik, igrd   ! dummy loop indices 
    7979      INTEGER  ::   ii, ij         ! 2D addresses 
     
    9292 
    9393 
    94    SUBROUTINE bdy_orl( idx, ptb, pta, dta, ll_npo ) 
     94   SUBROUTINE bdy_orl( idx, ptb, pta, dta, lrim0, ll_npo ) 
    9595      !!---------------------------------------------------------------------- 
    9696      !!                 ***  SUBROUTINE bdy_orl  *** 
     
    104104      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptb  ! before tracer field 
    105105      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta  ! tracer trend 
     106      LOGICAL                 , OPTIONAL,  INTENT(in) ::   lrim0   ! indicate if rim 0 is treated 
    106107      LOGICAL,                             INTENT(in) ::   ll_npo  ! switch for NPO version 
    107108      !! 
     
    111112      igrd = 1                       ! Everything is at T-points here 
    112113      ! 
    113       CALL bdy_orlanski_3d( idx, igrd, ptb(:,:,:), pta(:,:,:), dta, ll_npo ) 
     114      CALL bdy_orlanski_3d( idx, igrd, ptb(:,:,:), pta(:,:,:), dta, lrim0, ll_npo ) 
    114115      ! 
    115116   END SUBROUTINE bdy_orl 
    116117 
    117118 
    118    SUBROUTINE bdy_orlanski_2d( idx, igrd, phib, phia, phi_ext, ll_npo ) 
     119   SUBROUTINE bdy_orlanski_2d( idx, igrd, phib, phia, phi_ext, lrim0, ll_npo ) 
    119120      !!---------------------------------------------------------------------- 
    120121      !!                 ***  SUBROUTINE bdy_orlanski_2d  *** 
     
    132133      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   phia     ! model after 2D field (to be updated) 
    133134      REAL(wp), DIMENSION(:)  , INTENT(in   ) ::   phi_ext  ! external forcing data 
     135      LOGICAL, OPTIONAL,        INTENT(in   ) ::   lrim0    ! indicate if rim 0 is treated 
    134136      LOGICAL ,                 INTENT(in   ) ::   ll_npo   ! switch for NPO version 
    135137      ! 
     
    140142      INTEGER  ::   ii_offset, ij_offset                   ! offsets for mask indices 
    141143      INTEGER  ::   flagu, flagv                           ! short cuts 
     144      INTEGER  ::   ibeg, iend                             ! length of rim to be treated (rim 0 or rim 1 or both) 
    142145      REAL(wp) ::   zmask_x, zmask_y1, zmask_y2 
    143146      REAL(wp) ::   zex1, zex2, zey, zey1, zey2 
     
    146149      REAL(wp) ::   zdy_1, zdy_2, zsign_ups 
    147150      REAL(wp), PARAMETER :: zepsilon = 1.e-30                 ! local small value 
    148       REAL(wp), POINTER, DIMENSION(:,:)          :: pmask      ! land/sea mask for field 
    149       REAL(wp), POINTER, DIMENSION(:,:)          :: pmask_xdif ! land/sea mask for x-derivatives 
    150       REAL(wp), POINTER, DIMENSION(:,:)          :: pmask_ydif ! land/sea mask for y-derivatives 
     151      REAL(wp), POINTER, DIMENSION(:,:)          :: zmask      ! land/sea mask for field 
     152      REAL(wp), POINTER, DIMENSION(:,:)          :: zmask_xdif ! land/sea mask for x-derivatives 
     153      REAL(wp), POINTER, DIMENSION(:,:)          :: zmask_ydif ! land/sea mask for y-derivatives 
    151154      REAL(wp), POINTER, DIMENSION(:,:)          :: pe_xdif    ! scale factors for x-derivatives 
    152155      REAL(wp), POINTER, DIMENSION(:,:)          :: pe_ydif    ! scale factors for y-derivatives 
     
    159162      SELECT CASE(igrd) 
    160163         CASE(1) 
    161             pmask      => tmask(:,:,1) 
    162             pmask_xdif => umask(:,:,1) 
    163             pmask_ydif => vmask(:,:,1) 
     164            zmask      => tmask(:,:,1) 
     165            zmask_xdif => umask(:,:,1) 
     166            zmask_ydif => vmask(:,:,1) 
    164167            pe_xdif    => e1u(:,:) 
    165168            pe_ydif    => e2v(:,:) 
     
    167170            ij_offset = 0 
    168171         CASE(2) 
    169             pmask      => umask(:,:,1) 
    170             pmask_xdif => tmask(:,:,1) 
    171             pmask_ydif => fmask(:,:,1) 
     172            zmask      => umask(:,:,1) 
     173            zmask_xdif => tmask(:,:,1) 
     174            zmask_ydif => fmask(:,:,1) 
    172175            pe_xdif    => e1t(:,:) 
    173176            pe_ydif    => e2f(:,:) 
     
    175178            ij_offset = 0 
    176179         CASE(3) 
    177             pmask      => vmask(:,:,1) 
    178             pmask_xdif => fmask(:,:,1) 
    179             pmask_ydif => tmask(:,:,1) 
     180            zmask      => vmask(:,:,1) 
     181            zmask_xdif => fmask(:,:,1) 
     182            zmask_ydif => tmask(:,:,1) 
    180183            pe_xdif    => e1f(:,:) 
    181184            pe_ydif    => e2t(:,:) 
     
    185188      END SELECT 
    186189      ! 
    187       DO jb = 1, idx%nblenrim(igrd) 
     190      IF( PRESENT(lrim0) ) THEN 
     191         IF( lrim0 ) THEN   ;   ibeg = 1                       ;   iend = idx%nblenrim0(igrd)   ! rim 0 
     192         ELSE               ;   ibeg = idx%nblenrim0(igrd)+1   ;   iend = idx%nblenrim(igrd)    ! rim 1 
     193         END IF 
     194      ELSE                  ;   ibeg = 1                       ;   iend = idx%nblenrim(igrd)    ! both 
     195      END IF 
     196      ! 
     197      DO jb = ibeg, iend 
    188198         ii  = idx%nbi(jb,igrd) 
    189199         ij  = idx%nbj(jb,igrd)  
     200         IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )   CYCLE 
    190201         flagu = int( idx%flagu(jb,igrd) ) 
    191202         flagv = int( idx%flagv(jb,igrd) ) 
     
    203214         ! 
    204215         ! Calculate scale factors for calculation of spatial derivatives.         
    205          zex1 = ( abs(iibm1-iibm2) * pe_xdif(iibm1+ii_offset,ijbm1          )         & 
    206         &       + abs(ijbm1-ijbm2) * pe_ydif(iibm1          ,ijbm1+ij_offset) )  
    207          zex2 = ( abs(iibm1-iibm2) * pe_xdif(iibm2+ii_offset,ijbm2          )         & 
    208         &       + abs(ijbm1-ijbm2) * pe_ydif(iibm2          ,ijbm2+ij_offset) )  
    209          zey1 = ( (iibm1-iibm1jm1) * pe_xdif(iibm1jm1+ii_offset,ijbm1jm1          )  &  
     216         zex1 = ( abs(iibm1-iibm2) * pe_xdif(iibm1   +ii_offset,ijbm1             )   & 
     217        &       + abs(ijbm1-ijbm2) * pe_ydif(iibm1             ,ijbm1   +ij_offset) )  
     218         zex2 = ( abs(iibm1-iibm2) * pe_xdif(iibm2   +ii_offset,ijbm2             )   & 
     219        &       + abs(ijbm1-ijbm2) * pe_ydif(iibm2             ,ijbm2   +ij_offset) )  
     220         zey1 = ( (iibm1-iibm1jm1) * pe_xdif(iibm1jm1+ii_offset,ijbm1jm1          )   &  
    210221        &      +  (ijbm1-ijbm1jm1) * pe_ydif(iibm1jm1          ,ijbm1jm1+ij_offset) )  
    211          zey2 = ( (iibm1jp1-iibm1) * pe_xdif(iibm1+ii_offset,ijbm1)                  & 
    212         &      +  (ijbm1jp1-ijbm1) * pe_ydif(iibm1          ,ijbm1+ij_offset) )  
     222         zey2 = ( (iibm1jp1-iibm1) * pe_xdif(iibm1   +ii_offset,ijbm1             )   & 
     223        &      +  (ijbm1jp1-ijbm1) * pe_ydif(iibm1             ,ijbm1   +ij_offset) )  
    213224         ! make sure scale factors are nonzero 
    214225         if( zey1 .lt. rsmall ) zey1 = zey2 
     
    217228         zey1 = max(zey1,rsmall); zey2 = max(zey2,rsmall);  
    218229         ! 
    219          ! Calculate masks for calculation of spatial derivatives.         
    220          zmask_x = ( abs(iibm1-iibm2) * pmask_xdif(iibm2+ii_offset,ijbm2          )         & 
    221         &          + abs(ijbm1-ijbm2) * pmask_ydif(iibm2          ,ijbm2+ij_offset) )  
    222          zmask_y1 = ( (iibm1-iibm1jm1) * pmask_xdif(iibm1jm1+ii_offset,ijbm1jm1          )  &  
    223         &          +  (ijbm1-ijbm1jm1) * pmask_ydif(iibm1jm1          ,ijbm1jm1+ij_offset) )  
    224          zmask_y2 = ( (iibm1jp1-iibm1) * pmask_xdif(iibm1+ii_offset,ijbm1)                  & 
    225         &          +  (ijbm1jp1-ijbm1) * pmask_ydif(iibm1          ,ijbm1+ij_offset) )  
     230         ! Calculate masks for calculation of spatial derivatives. 
     231         zmask_x  = ( abs(iibm1-iibm2) * zmask_xdif(iibm2   +ii_offset,ijbm2               )   & 
     232        &           + abs(ijbm1-ijbm2) * zmask_ydif(iibm2             ,ijbm2   +ij_offset) )  
     233         zmask_y1 = ( (iibm1-iibm1jm1) * zmask_xdif(iibm1jm1+ii_offset,ijbm1jm1            )   &  
     234        &          +  (ijbm1-ijbm1jm1) * zmask_ydif(iibm1jm1          ,ijbm1jm1+ij_offset) )  
     235         zmask_y2 = ( (iibm1jp1-iibm1) * zmask_xdif(iibm1   +ii_offset,ijbm1               )   & 
     236        &          +  (ijbm1jp1-ijbm1) * zmask_ydif(iibm1             ,ijbm1   +ij_offset) )  
    226237 
    227238         ! Calculation of terms required for both versions of the scheme.  
     
    231242         ! Note no rdt factor in expression for zdt because it cancels in the expressions for  
    232243         ! zrx and zry. 
    233          zdt = phia(iibm1,ijbm1) - phib(iibm1,ijbm1) 
    234          zdx = ( ( phia(iibm1,ijbm1) - phia(iibm2,ijbm2) ) / zex2 ) * zmask_x  
     244         zdt   =     phia(iibm1   ,ijbm1   ) - phib(iibm1   ,ijbm1   ) 
     245         zdx   = ( ( phia(iibm1   ,ijbm1   ) - phia(iibm2   ,ijbm2   ) ) / zex2 ) * zmask_x  
    235246         zdy_1 = ( ( phib(iibm1   ,ijbm1   ) - phib(iibm1jm1,ijbm1jm1) ) / zey1 ) * zmask_y1     
    236          zdy_2 = ( ( phib(iibm1jp1,ijbm1jp1) - phib(iibm1   ,ijbm1)    ) / zey2 ) * zmask_y2  
     247         zdy_2 = ( ( phib(iibm1jp1,ijbm1jp1) - phib(iibm1   ,ijbm1   ) ) / zey2 ) * zmask_y2  
    237248         zdy_centred = 0.5 * ( zdy_1 + zdy_2 ) 
    238249!!$         zdy_centred = phib(iibm1jp1,ijbm1jp1) - phib(iibm1jm1,ijbm1jm1) 
     
    265276           &                    + zwgt * ( phi_ext(jb) - phib(ii,ij) ) ) / ( 1. + zrx )  
    266277         end if 
    267          phia(ii,ij) = phia(ii,ij) * pmask(ii,ij) 
     278         phia(ii,ij) = phia(ii,ij) * zmask(ii,ij) 
    268279      END DO 
    269280      ! 
     
    271282 
    272283 
    273    SUBROUTINE bdy_orlanski_3d( idx, igrd, phib, phia, phi_ext, ll_npo ) 
     284   SUBROUTINE bdy_orlanski_3d( idx, igrd, phib, phia, phi_ext, lrim0, ll_npo ) 
    274285      !!---------------------------------------------------------------------- 
    275286      !!                 ***  SUBROUTINE bdy_orlanski_3d  *** 
     
    287298      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phia     ! model after 3D field (to be updated) 
    288299      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   phi_ext  ! external forcing data 
     300      LOGICAL, OPTIONAL,          INTENT(in   ) ::   lrim0    ! indicate if rim 0 is treated 
    289301      LOGICAL ,                   INTENT(in   ) ::   ll_npo   ! switch for NPO version 
    290302      ! 
     
    295307      INTEGER  ::   ii_offset, ij_offset                   ! offsets for mask indices 
    296308      INTEGER  ::   flagu, flagv                           ! short cuts 
     309      INTEGER  ::   ibeg, iend                             ! length of rim to be treated (rim 0 or rim 1 or both) 
    297310      REAL(wp) ::   zmask_x, zmask_y1, zmask_y2 
    298311      REAL(wp) ::   zex1, zex2, zey, zey1, zey2 
     
    301314      REAL(wp) ::   zdy_1, zdy_2,  zsign_ups 
    302315      REAL(wp), PARAMETER :: zepsilon = 1.e-30                 ! local small value 
    303       REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask      ! land/sea mask for field 
    304       REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask_xdif ! land/sea mask for x-derivatives 
    305       REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask_ydif ! land/sea mask for y-derivatives 
     316      REAL(wp), POINTER, DIMENSION(:,:,:)        :: zmask      ! land/sea mask for field 
     317      REAL(wp), POINTER, DIMENSION(:,:,:)        :: zmask_xdif ! land/sea mask for x-derivatives 
     318      REAL(wp), POINTER, DIMENSION(:,:,:)        :: zmask_ydif ! land/sea mask for y-derivatives 
    306319      REAL(wp), POINTER, DIMENSION(:,:)          :: pe_xdif    ! scale factors for x-derivatives 
    307320      REAL(wp), POINTER, DIMENSION(:,:)          :: pe_ydif    ! scale factors for y-derivatives 
     
    314327      SELECT CASE(igrd) 
    315328         CASE(1) 
    316             pmask      => tmask(:,:,:) 
    317             pmask_xdif => umask(:,:,:) 
    318             pmask_ydif => vmask(:,:,:) 
     329            zmask      => tmask(:,:,:) 
     330            zmask_xdif => umask(:,:,:) 
     331            zmask_ydif => vmask(:,:,:) 
    319332            pe_xdif    => e1u(:,:) 
    320333            pe_ydif    => e2v(:,:) 
     
    322335            ij_offset = 0 
    323336         CASE(2) 
    324             pmask      => umask(:,:,:) 
    325             pmask_xdif => tmask(:,:,:) 
    326             pmask_ydif => fmask(:,:,:) 
     337            zmask      => umask(:,:,:) 
     338            zmask_xdif => tmask(:,:,:) 
     339            zmask_ydif => fmask(:,:,:) 
    327340            pe_xdif    => e1t(:,:) 
    328341            pe_ydif    => e2f(:,:) 
     
    330343            ij_offset = 0 
    331344         CASE(3) 
    332             pmask      => vmask(:,:,:) 
    333             pmask_xdif => fmask(:,:,:) 
    334             pmask_ydif => tmask(:,:,:) 
     345            zmask      => vmask(:,:,:) 
     346            zmask_xdif => fmask(:,:,:) 
     347            zmask_ydif => tmask(:,:,:) 
    335348            pe_xdif    => e1f(:,:) 
    336349            pe_ydif    => e2t(:,:) 
     
    339352         CASE DEFAULT ;   CALL ctl_stop( 'unrecognised value for igrd in bdy_orlanksi_2d' ) 
    340353      END SELECT 
    341  
     354      ! 
     355      IF( PRESENT(lrim0) ) THEN 
     356         IF( lrim0 ) THEN   ;   ibeg = 1                       ;   iend = idx%nblenrim0(igrd)   ! rim 0 
     357         ELSE               ;   ibeg = idx%nblenrim0(igrd)+1   ;   iend = idx%nblenrim(igrd)    ! rim 1 
     358         END IF 
     359      ELSE                  ;   ibeg = 1                       ;   iend = idx%nblenrim(igrd)    ! both 
     360      END IF 
     361      ! 
    342362      DO jk = 1, jpk 
    343363         !             
    344          DO jb = 1, idx%nblenrim(igrd) 
     364         DO jb = ibeg, iend 
    345365            ii  = idx%nbi(jb,igrd) 
    346366            ij  = idx%nbj(jb,igrd)  
     367            IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )   CYCLE 
    347368            flagu = int( idx%flagu(jb,igrd) ) 
    348369            flagv = int( idx%flagv(jb,igrd) ) 
     
    360381            ! 
    361382            ! Calculate scale factors for calculation of spatial derivatives.         
    362             zex1 = ( abs(iibm1-iibm2) * pe_xdif(iibm1+ii_offset,ijbm1          )         & 
    363            &       + abs(ijbm1-ijbm2) * pe_ydif(iibm1          ,ijbm1+ij_offset) )  
    364             zex2 = ( abs(iibm1-iibm2) * pe_xdif(iibm2+ii_offset,ijbm2          )         & 
    365            &       + abs(ijbm1-ijbm2) * pe_ydif(iibm2          ,ijbm2+ij_offset) )  
    366             zey1 = ( (iibm1-iibm1jm1) * pe_xdif(iibm1jm1+ii_offset,ijbm1jm1          )  &  
     383            zex1 = ( abs(iibm1-iibm2) * pe_xdif(iibm1   +ii_offset,ijbm1             )   & 
     384           &       + abs(ijbm1-ijbm2) * pe_ydif(iibm1             ,ijbm1+ij_offset   ) )  
     385            zex2 = ( abs(iibm1-iibm2) * pe_xdif(iibm2   +ii_offset,ijbm2             )   & 
     386           &       + abs(ijbm1-ijbm2) * pe_ydif(iibm2             ,ijbm2+ij_offset   ) )  
     387            zey1 = ( (iibm1-iibm1jm1) * pe_xdif(iibm1jm1+ii_offset,ijbm1jm1          )   &  
    367388           &      +  (ijbm1-ijbm1jm1) * pe_ydif(iibm1jm1          ,ijbm1jm1+ij_offset) )  
    368             zey2 = ( (iibm1jp1-iibm1) * pe_xdif(iibm1+ii_offset,ijbm1)                  & 
    369            &      +  (ijbm1jp1-ijbm1) * pe_ydif(iibm1          ,ijbm1+ij_offset) )  
     389            zey2 = ( (iibm1jp1-iibm1) * pe_xdif(iibm1   +ii_offset,ijbm1             )   & 
     390           &      +  (ijbm1jp1-ijbm1) * pe_ydif(iibm1             ,ijbm1+ij_offset   ) )  
    370391            ! make sure scale factors are nonzero 
    371392            if( zey1 .lt. rsmall ) zey1 = zey2 
     
    375396            ! 
    376397            ! Calculate masks for calculation of spatial derivatives.         
    377             zmask_x = ( abs(iibm1-iibm2) * pmask_xdif(iibm2+ii_offset,ijbm2          ,jk)          & 
    378            &          + abs(ijbm1-ijbm2) * pmask_ydif(iibm2          ,ijbm2+ij_offset,jk) )  
    379             zmask_y1 = ( (iibm1-iibm1jm1) * pmask_xdif(iibm1jm1+ii_offset,ijbm1jm1          ,jk)   &  
    380            &          +  (ijbm1-ijbm1jm1) * pmask_ydif(iibm1jm1          ,ijbm1jm1+ij_offset,jk) )  
    381             zmask_y2 = ( (iibm1jp1-iibm1) * pmask_xdif(iibm1+ii_offset,ijbm1          ,jk)         & 
    382            &          +  (ijbm1jp1-ijbm1) * pmask_ydif(iibm1          ,ijbm1+ij_offset,jk) )  
     398            zmask_x  = ( abs(iibm1-iibm2) * zmask_xdif(iibm2   +ii_offset,ijbm2             ,jk)   & 
     399           &           + abs(ijbm1-ijbm2) * zmask_ydif(iibm2             ,ijbm2   +ij_offset,jk) )  
     400            zmask_y1 = ( (iibm1-iibm1jm1) * zmask_xdif(iibm1jm1+ii_offset,ijbm1jm1          ,jk)   &  
     401           &          +  (ijbm1-ijbm1jm1) * zmask_ydif(iibm1jm1          ,ijbm1jm1+ij_offset,jk) )  
     402            zmask_y2 = ( (iibm1jp1-iibm1) * zmask_xdif(iibm1   +ii_offset,ijbm1             ,jk)   & 
     403           &          +  (ijbm1jp1-ijbm1) * zmask_ydif(iibm1             ,ijbm1   +ij_offset,jk) )  
    383404            ! 
    384405            ! Calculate normal (zrx) and tangential (zry) components of radiation velocities. 
     
    386407            ! Centred derivative is calculated as average of "left" and "right" derivatives for  
    387408            ! this reason.  
    388             zdt = phia(iibm1,ijbm1,jk) - phib(iibm1,ijbm1,jk) 
    389             zdx = ( ( phia(iibm1,ijbm1,jk) - phia(iibm2,ijbm2,jk) ) / zex2 ) * zmask_x                   
     409            zdt   =     phia(iibm1   ,ijbm1   ,jk) - phib(iibm1   ,ijbm1   ,jk) 
     410            zdx   = ( ( phia(iibm1   ,ijbm1   ,jk) - phia(iibm2   ,ijbm2   ,jk) ) / zex2 ) * zmask_x                   
    390411            zdy_1 = ( ( phib(iibm1   ,ijbm1   ,jk) - phib(iibm1jm1,ijbm1jm1,jk) ) / zey1 ) * zmask_y1   
    391412            zdy_2 = ( ( phib(iibm1jp1,ijbm1jp1,jk) - phib(iibm1   ,ijbm1   ,jk) ) / zey2 ) * zmask_y2       
     
    421442              &                       + zwgt * ( phi_ext(jb,jk) - phib(ii,ij,jk) ) ) / ( 1. + zrx )  
    422443            end if 
    423             phia(ii,ij,jk) = phia(ii,ij,jk) * pmask(ii,ij,jk) 
     444            phia(ii,ij,jk) = phia(ii,ij,jk) * zmask(ii,ij,jk) 
    424445         END DO 
    425446         ! 
     
    428449   END SUBROUTINE bdy_orlanski_3d 
    429450 
    430    SUBROUTINE bdy_nmn( idx, igrd, phia ) 
     451   SUBROUTINE bdy_nmn( idx, igrd, phia, lrim0 ) 
    431452      !!---------------------------------------------------------------------- 
    432453      !!                 ***  SUBROUTINE bdy_nmn  *** 
     
    434455      !! ** Purpose : Duplicate the value at open boundaries, zero gradient. 
    435456      !!  
    436       !!---------------------------------------------------------------------- 
    437       INTEGER,                    INTENT(in)     ::   igrd     ! grid index 
    438       REAL(wp), DIMENSION(:,:,:), INTENT(inout)  ::   phia     ! model after 3D field (to be updated) 
    439       TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices 
     457      !! 
     458      !! ** Method  : - take the average of free ocean neighbours 
     459      !! 
     460      !!      ___   !   |_____|   !   ___|    !   __|x o   !   |_   _|     ! |       
     461      !!   __|x     !      x      !     x o   !      o     !     |_|       ! |x o    
     462      !!      o     !      o      !     o     !            !    o x o      ! |x_x_  
     463      !!                                                   !      o       
     464      !!---------------------------------------------------------------------- 
     465      INTEGER,                    INTENT(in   )  ::   igrd     ! grid index 
     466      REAL(wp), DIMENSION(:,:,:), INTENT(inout)  ::   phia     ! model after 3D field (to be updated), must be masked 
     467      TYPE(OBC_INDEX),            INTENT(in   )  ::   idx      ! OBC indices 
     468      LOGICAL, OPTIONAL,          INTENT(in   )  ::   lrim0    ! indicate if rim 0 is treated 
    440469      !!  
    441       REAL(wp) ::   zcoef, zcoef1, zcoef2 
    442       REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask      ! land/sea mask for field 
    443       REAL(wp), POINTER, DIMENSION(:,:)        :: bdypmask      ! land/sea mask for field 
     470      REAL(wp) ::   zweight 
     471      REAL(wp), POINTER, DIMENSION(:,:,:)      :: zmask         ! land/sea mask for field 
    444472      INTEGER  ::   ib, ik   ! dummy loop indices 
    445       INTEGER  ::   ii, ij, ip, jp   ! 2D addresses 
    446       !!---------------------------------------------------------------------- 
     473      INTEGER  ::   ii, ij   ! 2D addresses 
     474      INTEGER  ::   ipkm1    ! size of phia third dimension minus 1 
     475      INTEGER  ::   ibeg, iend                          ! length of rim to be treated (rim 0 or rim 1 or both) 
     476      INTEGER  ::   ii1, ii2, ii3, ij1, ij2, ij3, itreat 
     477      !!---------------------------------------------------------------------- 
     478      ! 
     479      ipkm1 = MAX( SIZE(phia,3) - 1, 1 )  
    447480      ! 
    448481      SELECT CASE(igrd) 
    449          CASE(1) 
    450             pmask => tmask(:,:,:) 
    451             bdypmask => bdytmask(:,:) 
    452          CASE(2) 
    453             pmask => umask(:,:,:) 
    454             bdypmask => bdyumask(:,:) 
    455          CASE(3) 
    456             pmask => vmask(:,:,:) 
    457             bdypmask => bdyvmask(:,:) 
     482         CASE(1)   ;   zmask => tmask(:,:,:) 
     483         CASE(2)   ;   zmask => umask(:,:,:) 
     484         CASE(3)   ;   zmask => vmask(:,:,:) 
    458485         CASE DEFAULT ;   CALL ctl_stop( 'unrecognised value for igrd in bdy_nmn' ) 
    459486      END SELECT 
    460       DO ib = 1, idx%nblenrim(igrd) 
     487      ! 
     488      IF( PRESENT(lrim0) ) THEN 
     489         IF( lrim0 ) THEN   ;   ibeg = 1                       ;   iend = idx%nblenrim0(igrd)   ! rim 0 
     490         ELSE               ;   ibeg = idx%nblenrim0(igrd)+1   ;   iend = idx%nblenrim(igrd)    ! rim 1 
     491         END IF 
     492      ELSE                  ;   ibeg = 1                       ;   iend = idx%nblenrim(igrd)    ! both 
     493      END IF 
     494      ! 
     495      DO ib = ibeg, iend 
    461496         ii = idx%nbi(ib,igrd) 
    462497         ij = idx%nbj(ib,igrd) 
    463          DO ik = 1, jpkm1 
    464             ! search the sense of the gradient 
    465             zcoef1 = bdypmask(ii-1,ij  )*pmask(ii-1,ij,ik) +  bdypmask(ii+1,ij  )*pmask(ii+1,ij,ik) 
    466             zcoef2 = bdypmask(ii  ,ij-1)*pmask(ii,ij-1,ik) +  bdypmask(ii  ,ij+1)*pmask(ii,ij+1,ik) 
    467             IF ( nint(zcoef1+zcoef2) == 0) THEN 
    468                ! corner **** we probably only want to set the tangentail component for the dynamics here 
    469                zcoef = pmask(ii-1,ij,ik) + pmask(ii+1,ij,ik) +  pmask(ii,ij-1,ik) +  pmask(ii,ij+1,ik) 
    470                IF (zcoef > .5_wp) THEN ! Only set none isolated points. 
    471                  phia(ii,ij,ik) = phia(ii-1,ij  ,ik) * pmask(ii-1,ij  ,ik) + & 
    472                    &              phia(ii+1,ij  ,ik) * pmask(ii+1,ij  ,ik) + & 
    473                    &              phia(ii  ,ij-1,ik) * pmask(ii  ,ij-1,ik) + & 
    474                    &              phia(ii  ,ij+1,ik) * pmask(ii  ,ij+1,ik) 
    475                  phia(ii,ij,ik) = ( phia(ii,ij,ik) / zcoef ) * pmask(ii,ij,ik) 
    476                ELSE 
    477                  phia(ii,ij,ik) = phia(ii,ij  ,ik) * pmask(ii,ij  ,ik) 
    478                ENDIF 
    479             ELSEIF ( nint(zcoef1+zcoef2) == 2) THEN 
    480                ! oblique corner **** we probably only want to set the normal component for the dynamics here 
    481                zcoef = pmask(ii-1,ij,ik)*bdypmask(ii-1,ij  ) + pmask(ii+1,ij,ik)*bdypmask(ii+1,ij  ) + & 
    482                    &   pmask(ii,ij-1,ik)*bdypmask(ii,ij -1 ) +  pmask(ii,ij+1,ik)*bdypmask(ii,ij+1  ) 
    483                phia(ii,ij,ik) = phia(ii-1,ij  ,ik) * pmask(ii-1,ij  ,ik)*bdypmask(ii-1,ij  ) + & 
    484                    &            phia(ii+1,ij  ,ik) * pmask(ii+1,ij  ,ik)*bdypmask(ii+1,ij  )  + & 
    485                    &            phia(ii  ,ij-1,ik) * pmask(ii  ,ij-1,ik)*bdypmask(ii,ij -1 ) + & 
    486                    &            phia(ii  ,ij+1,ik) * pmask(ii  ,ij+1,ik)*bdypmask(ii,ij+1  ) 
    487   
    488                phia(ii,ij,ik) = ( phia(ii,ij,ik) / MAX(1._wp, zcoef) ) * pmask(ii,ij,ik) 
    489             ELSE 
    490                ip = nint(bdypmask(ii+1,ij  )*pmask(ii+1,ij,ik) - bdypmask(ii-1,ij  )*pmask(ii-1,ij,ik)) 
    491                jp = nint(bdypmask(ii  ,ij+1)*pmask(ii,ij+1,ik) - bdypmask(ii  ,ij-1)*pmask(ii,ij-1,ik)) 
    492                phia(ii,ij,ik) = phia(ii+ip,ij+jp,ik) * pmask(ii+ip,ij+jp,ik) 
    493             ENDIF 
    494          END DO 
     498         itreat = idx%ntreat(ib,igrd) 
     499         CALL find_neib( ii, ij, itreat, ii1, ij1, ii2, ij2, ii3, ij3 )   ! find free ocean neighbours 
     500         SELECT CASE( itreat ) 
     501         CASE( 1:8 ) 
     502            IF( ii1 < 1 .OR. ii1 > jpi .OR. ij1 < 1 .OR. ij1 > jpj )   CYCLE 
     503            DO ik = 1, ipkm1 
     504               IF( zmask(ii1,ij1,ik) /= 0. )   phia(ii,ij,ik) = phia(ii1,ij1,ik)   
     505            END DO 
     506         CASE( 9:12 ) 
     507            IF( ii1 < 1 .OR. ii1 > jpi .OR. ij1 < 1 .OR. ij1 > jpj )   CYCLE 
     508            IF( ii2 < 1 .OR. ii2 > jpi .OR. ij2 < 1 .OR. ij2 > jpj )   CYCLE 
     509            DO ik = 1, ipkm1 
     510               zweight = zmask(ii1,ij1,ik) + zmask(ii2,ij2,ik) 
     511               IF( zweight /= 0. )   phia(ii,ij,ik) = ( phia(ii1,ij1,ik) + phia(ii2,ij2,ik) ) / zweight 
     512            END DO 
     513         CASE( 13:16 ) 
     514            IF( ii1 < 1 .OR. ii1 > jpi .OR. ij1 < 1 .OR. ij1 > jpj )   CYCLE 
     515            IF( ii2 < 1 .OR. ii2 > jpi .OR. ij2 < 1 .OR. ij2 > jpj )   CYCLE 
     516            IF( ii3 < 1 .OR. ii3 > jpi .OR. ij3 < 1 .OR. ij3 > jpj )   CYCLE 
     517            DO ik = 1, ipkm1 
     518               zweight = zmask(ii1,ij1,ik) + zmask(ii2,ij2,ik) + zmask(ii3,ij3,ik) 
     519               IF( zweight /= 0. )   phia(ii,ij,ik) = ( phia(ii1,ij1,ik) + phia(ii2,ij2,ik) + phia(ii3,ij3,ik) ) / zweight 
     520            END DO 
     521         END SELECT 
    495522      END DO 
    496523      ! 
Note: See TracChangeset for help on using the changeset viewer.