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 11183 for NEMO – NEMO

Changeset 11183 for NEMO


Ignore:
Timestamp:
2019-06-26T16:09:28+02:00 (5 years ago)
Author:
girrmann
Message:

dev_r10984_HPC-13 : improvement of neumann boundary condition, changes the results, see #2285

Location:
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyini.F90

    r11071 r11183  
    3434 
    3535   PUBLIC   bdy_init   ! routine called in nemo_init 
     36   PUBLIC   find_neib   ! routine called in bdy_nmn 
    3637 
    3738   INTEGER, PARAMETER ::   jp_nseg = 100   !  
     
    131132      INTEGER  ::   jpbdtau, jpbdtas                       !   -       - 
    132133      INTEGER  ::   ib_bdy1, ib_bdy2, ib1, ib2             !   -       - 
    133       INTEGER  ::   i_offset, j_offset, inbdy, itreat      !   -       - 
    134       INTEGER  ::   ii1, ii2, ii3, ij1, ij2, ij3, iibe, ijbe   !   -       - 
    135       INTEGER  ::   flagu, flagv                               ! short cuts 
     134      INTEGER  ::   i_offset, j_offset, inn                !   -       - 
     135      INTEGER  ::   ii1, ii2, ii3, ij1, ij2, ij3           !   -       - 
     136      INTEGER  ::   iibe, ijbe, iibi, ijbi                 !   -       - 
     137      INTEGER  ::   flagu, flagv                           ! short cuts 
    136138      INTEGER , POINTER  ::  nbi, nbj, nbr                 ! short cuts 
    137       REAL(wp), POINTER, DIMENSION(:,:)       ::   pmask    ! pointer to 2D mask fields 
     139      REAL(wp), POINTER, DIMENSION(:,:)       ::   pmask   ! pointer to 2D mask fields 
    138140      REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars 
    139141      INTEGER, DIMENSION (2)                  ::   kdimsz 
     
    147149      REAL(wp), TARGET, DIMENSION(jpi,jpj) ::   zfmask  ! temporary fmask array excluding coastal boundary condition (shlat) 
    148150      REAL(wp)        , DIMENSION(jpi,jpj) ::   ztmp 
    149       LOGICAL  ::   llnobdy, llsobdy, lleabdy, llwebdy     ! local logicals 
     151      LOGICAL  ::   llnon, llson, llean, llwen          ! local logicals 
    150152      !! 
    151153      CHARACTER(LEN=1)                     ::   ctypebdy   !     -        -  
     
    10721074      END DO 
    10731075       
    1074       ! detect corner interior and its orientation index 1 to 4  depending on the orientation 
    1075       ! detect corner exterior and its orientation index 5 to 8  depending on the orientation 
    1076       ! detect geometries with 3 neighbours        index 9 to 12 depending on the orientation 
    1077       ! else                                       index 0 
    10781076      DO ib_bdy = 1, nb_bdy 
    10791077         DO igrd = 1, jpbgrd 
     
    10871085               ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    10881086               ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    1089                IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )   CYCLE 
    1090                llnobdy = pmask(ii  ,ij+1) == 1.   
    1091                llsobdy = pmask(ii  ,ij-1) == 1.  
    1092                lleabdy = pmask(ii+1,ij  ) == 1.  
    1093                llwebdy = pmask(ii-1,ij  ) == 1.  
    1094                inbdy  = COUNT( (/ llnobdy, llsobdy, lleabdy, llwebdy /) ) 
    1095                IF( inbdy == 0 )   THEN   ! no neighbours -> interior of a corner 
    1096                   !               !              !     _____     !     _____      
    1097                   !  1 |   o      !  2  o   |    !  3 | x        !  4     x |     
    1098                   !    |_x_ _     !    _ _x_|    !    |   o      !      o   | 
    1099                   IF( pmask(ii+1,ij+1) == 1. )   ztmp(ii,ij) = 1. 
    1100                   IF( pmask(ii-1,ij+1) == 1. )   ztmp(ii,ij) = 2. 
    1101                   IF( pmask(ii+1,ij-1) == 1. )   ztmp(ii,ij) = 3. 
    1102                   IF( pmask(ii-1,ij-1) == 1. )   ztmp(ii,ij) = 4. 
     1087               IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )   CYCLE  
     1088               llean = pmask(ii+1,ij  ) == 1. 
     1089               llwen = pmask(ii-1,ij  ) == 1. 
     1090               llnon = pmask(ii  ,ij+1) == 1. 
     1091               llson = pmask(ii  ,ij-1) == 1. 
     1092               inn  = COUNT( (/ llnon, llson, llean, llwen /) ) 
     1093               IF( inn == 0 ) THEN   ! no neighbours -> interior of a corner  or  cluster of rim points 
     1094                  !               !              !     _____     !     _____    !    __     __ 
     1095                  !  1 |   o      !  2  o   |    !  3 | x        !  4     x |   !      |   |   -> error 
     1096                  !    |_x_ _     !    _ _x_|    !    |   o      !      o   |   !      |x_x| 
     1097                  IF(     pmask(ii+1,ij+1) == 1. ) THEN   ;   ztmp(ii,ij) = 1. 
     1098                  ELSEIF( pmask(ii-1,ij+1) == 1. ) THEN   ;   ztmp(ii,ij) = 2. 
     1099                  ELSEIF( pmask(ii+1,ij-1) == 1. ) THEN   ;   ztmp(ii,ij) = 3. 
     1100                  ELSEIF( pmask(ii-1,ij-1) == 1. ) THEN   ;   ztmp(ii,ij) = 4. 
     1101                  ELSE 
     1102                     WRITE(ctmp1,*) ' E R R O R : Problem with  ',cgrid(igrd) ,' grid point', ii, ij,   & 
     1103                          ' on boundary set ', ib_bdy, ' has no free ocean neighbour' 
     1104                     WRITE(ctmp2,*) ' There seems to be a cluster of rim points.' 
     1105                     WRITE(ctmp3,*) ' ========== ' 
     1106                     CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ' ) 
     1107                  END IF 
    11031108               END IF 
    1104                IF( inbdy == 1 )   THEN   ! middle of linear bdy 
    1105                   ztmp(ii,ij) = 0.   ! regular treatment with flags 
     1109               IF( inn == 1 ) THEN   ! middle of linear bdy  or incomplete corner  ! ___ o 
     1110                  !    |         !         |   !      o     !    ______            !    |x___ 
     1111                  ! 5  | x o     ! 6   o x |   ! 7  __x__   ! 8    x 
     1112                  !    |         !         |   !            !      o 
     1113                  IF( llean )   ztmp(ii,ij) = 5. 
     1114                  IF( llwen )   ztmp(ii,ij) = 6. 
     1115                  IF( llnon )   ztmp(ii,ij) = 7. 
     1116                  IF( llson )   ztmp(ii,ij) = 8. 
    11061117               END IF 
    1107                IF( inbdy == 2 )  THEN   ! exterior of a corner 
     1118               IF( inn == 2 ) THEN   ! exterior of a corner 
    11081119                  !        o      !        o      !    _____|       !       |_____   
    1109                   !  5 ____x o    !  6   o x___   ! 7      x o      !  8   o x       
     1120                  !  9 ____x o    ! 10   o x___   ! 11     x o      ! 12   o x       
    11101121                  !         |     !       |       !        o        !        o  
    1111                   IF( llnobdy .AND. lleabdy )   ztmp(ii,ij) = 5. 
    1112                   IF( llnobdy .AND. llwebdy )   ztmp(ii,ij) = 6. 
    1113                   IF( llsobdy .AND. lleabdy )   ztmp(ii,ij) = 7. 
    1114                   IF( llsobdy .AND. llwebdy )   ztmp(ii,ij) = 8. 
     1122                  IF( llnon .AND. llean )   ztmp(ii,ij) =  9. 
     1123                  IF( llnon .AND. llwen )   ztmp(ii,ij) = 10. 
     1124                  IF( llson .AND. llean )   ztmp(ii,ij) = 11. 
     1125                  IF( llson .AND. llwen )   ztmp(ii,ij) = 12. 
    11151126               END IF 
    1116                IF( inbdy == 3 )   THEN   ! 3 neighbours __   __ 
     1127               IF( inn == 3 ) THEN   ! 3 neighbours    __   __ 
    11171128                  !    |_  o      !        o  _|  !       |_|     !       o          
    1118                   !  9  _| x o    ! 10   o x |_   ! 11   o x o    ! 12  o x o        
     1129                  ! 13  _| x o    ! 14   o x |_   ! 15   o x o    ! 16  o x o        
    11191130                  !    |   o      !        o   |  !        o      !    __|¨|__     
    1120                   IF( llnobdy .AND. lleabdy .AND. llsobdy )   ztmp(ii,ij) = 9. 
    1121                   IF( llnobdy .AND. llwebdy .AND. llsobdy )   ztmp(ii,ij) = 10. 
    1122                   IF( llwebdy .AND. llsobdy .AND. lleabdy )   ztmp(ii,ij) = 11. 
    1123                   IF( llwebdy .AND. llnobdy .AND. lleabdy )   ztmp(ii,ij) = 12. 
     1131                  IF( llnon .AND. llean .AND. llson )   ztmp(ii,ij) = 13. 
     1132                  IF( llnon .AND. llwen .AND. llson )   ztmp(ii,ij) = 14. 
     1133                  IF( llwen .AND. llson .AND. llean )   ztmp(ii,ij) = 15. 
     1134                  IF( llwen .AND. llnon .AND. llean )   ztmp(ii,ij) = 16. 
    11241135               END IF 
    1125                IF( inbdy == 4 )  THEN 
    1126                   WRITE(ctmp1,*) ' E R R O R : Problem with  ',cgrid(igrd) ,' grid points,',   & 
    1127                        '  some points on boundary set ', ib_bdy, ' have 4 neighbours' 
     1136               IF( inn == 4 ) THEN 
     1137                  WRITE(ctmp1,*)  ' E R R O R : Problem with  ',cgrid(igrd) ,' grid point', ii, ij,   & 
     1138                       ' on boundary set ', ib_bdy, ' have 4 neighbours' 
    11281139                  WRITE(ctmp2,*) ' ========== ' 
    11291140                  CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     
    11311142            END DO 
    11321143            CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 
    1133             DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   
     1144            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
    11341145               ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    11351146               ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     
    11561167               flagv = NINT(idx_bdy(ib_bdy)%flagv(ib,igrd)) 
    11571168               iibe = ii - flagu   ! neighbouring point towards the exterior of the computational domain 
    1158                ijbe = ij - flagv   ! neighbouring point towards the exterior of the computational domain 
    1159                SELECT CASE( idx_bdy(ib_bdy)%ntreat(ib,igrd) )   ! points that will be used by bdy routines, -1 will be discarded 
    1160                CASE( 0 )    ;   ii1 = ii + flagu   ;   ii2 = -1   ;   ii3 = -1 
    1161                                 ij1 = ij + flagv   ;   ij2 = -1   ;   ij3 = -1 
    1162                CASE( 1  )   ;   ii1 = ii+1   ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1 
    1163                CASE( 2  )   ;   ii1 = ii-1   ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1 
    1164                CASE( 3  )   ;   ii1 = ii+1   ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1 
    1165                CASE( 4  )   ;   ii1 = ii-1   ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1 
    1166                CASE( 5  )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1 
    1167                CASE( 6  )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1 
    1168                CASE( 7  )   ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1  
    1169                CASE( 8  )   ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1 
    1170                CASE( 9  )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = ii     ;   ij3 = ij-1    
    1171                CASE( 10 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = ii     ;   ij3 = ij-1  
    1172                CASE( 11 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij-1   ;   ii3 = ii+1   ;   ij3 = ij    
    1173                CASE( 12 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij+1   ;   ii3 = ii+1   ;   ij3 = ij 
    1174                END SELECT 
     1169               ijbe = ij - flagv 
     1170               iibi = ii + flagu   ! neighbouring point towards the interior of the computational domain 
     1171               ijbi = ij + flagv 
     1172               CALL find_neib( ii, ij, idx_bdy(ib_bdy)%ntreat(ib,igrd), ii1, ij1, ii2, ij2, ii3, ij3 ) 
    11751173               ! 
    11761174               ! search neighbour in the  west/east  direction 
     
    11791177               ! (1)  o|x         OR    (2)   x|o 
    11801178               !       |___                 ___|  
    1181                IF( ii1 == 0     .OR. ii2 == 0     .OR. ii3 == 0     )   lrecv_bdyint(ib_bdy,igrd,1) = .true. 
    1182                IF( ii1 == jpi+1 .OR. ii2 == jpi+1 .OR. ii3 == jpi+1 )   lrecv_bdyint(ib_bdy,igrd,2) = .true.   
     1179               IF( iibi == 0     .OR. ii1 == 0     .OR. ii2 == 0     .OR. ii3 == 0     )   lrecv_bdyint(ib_bdy,igrd,1) = .true. 
     1180               IF( iibi == jpi+1 .OR. ii1 == jpi+1 .OR. ii2 == jpi+1 .OR. ii3 == jpi+1 )   lrecv_bdyint(ib_bdy,igrd,2) = .true.   
    11831181               IF( iibe == 0     )   lrecv_bdyext(ib_bdy,igrd,1) = .true. 
    11841182               IF( iibe == jpi+1 )   lrecv_bdyext(ib_bdy,igrd,2) = .true.   
     
    11871185               ! :     |  x:o    |    neighbour limited by ... would need o    |    o:x  |     : 
    11881186               ! :.....|_._:_____|   (1) W neighbour         E neighbour (2)   |_____:_._|.....: 
    1189                IF( ii == 2     .AND. ( nbondi ==  1 .OR. nbondi == 0 ) & 
    1190                     &          .AND. ( ii1 == 3     .OR. ii2 == 3     .OR. ii3 == 3     ) )   lsend_bdyint(ib_bdy,igrd,1) = .true. 
    1191                IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) & 
    1192                     &          .AND. ( ii1 == jpi-2 .OR. ii2 == jpi-2 .OR. ii3 == jpi-2 ) )   lsend_bdyint(ib_bdy,igrd,2) = .true. 
     1187               IF( ii == 2     .AND. ( nbondi ==  1 .OR. nbondi == 0 ) .AND. & 
     1188                    & ( iibi == 3     .OR. ii1 == 3     .OR. ii2 == 3     .OR. ii3 == 3    ) )   lsend_bdyint(ib_bdy,igrd,1) = .true. 
     1189               IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. & 
     1190                    & ( iibi == jpi-2 .OR. ii1 == jpi-2 .OR. ii2 == jpi-2 .OR. ii3 == jpi-2) )   lsend_bdyint(ib_bdy,igrd,2) = .true. 
    11931191               IF( ii == 2     .AND. ( nbondi ==  1 .OR. nbondi == 0 ) .AND. iibe == 3     )   lsend_bdyext(ib_bdy,igrd,1) = .true. 
    11941192               IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. iibe == jpi-2 )   lsend_bdyext(ib_bdy,igrd,2) = .true. 
     
    11991197               !  |   |___x___|   OR    |  |   x   | 
    12001198               !  v       o           (4)  |       | 
    1201                IF( ij1 == 0     .OR. ij2 == 0     .OR. ij3 == 0     )   lrecv_bdyint(ib_bdy,igrd,3) = .true. 
    1202                IF( ij1 == jpj+1 .OR. ij2 == jpj+1 .OR. ij3 == jpj+1 )   lrecv_bdyint(ib_bdy,igrd,4) = .true. 
     1199               IF( ijbi == 0     .OR. ij1 == 0     .OR. ij2 == 0     .OR. ij3 == 0     )   lrecv_bdyint(ib_bdy,igrd,3) = .true. 
     1200               IF( ijbi == jpj+1 .OR. ij1 == jpj+1 .OR. ij2 == jpj+1 .OR. ij3 == jpj+1 )   lrecv_bdyint(ib_bdy,igrd,4) = .true. 
    12031201               IF( ijbe == 0     )   lrecv_bdyext(ib_bdy,igrd,3) = .true. 
    12041202               IF( ijbe == jpj+1 )   lrecv_bdyext(ib_bdy,igrd,4) = .true. 
     
    12071205               !   |  |¨¨¨¨x¨¨¨¨|   neighbour limited by ... would need o     |  |....x....| 
    12081206               !      :_________:  (3) S neighbour          N neighbour (4)   v  |    o    |    
    1209                IF( ij == 2     .AND. ( nbondj ==  1 .OR. nbondj == 0 ) & 
    1210                     &          .AND. ( ij1 == 3     .OR. ij2 == 3     .OR. ij3 == 3     ) )   lsend_bdyint(ib_bdy,igrd,3) = .true. 
    1211                IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) & 
    1212                     &          .AND. ( ij1 == jpj-2 .OR. ij2 == jpj-2 .OR. ij3 == jpj-2 ) )   lsend_bdyint(ib_bdy,igrd,4) = .true. 
     1207               IF( ij == 2     .AND. ( nbondj ==  1 .OR. nbondj == 0 ) .AND. & 
     1208                    & ( ijbi == 3     .OR. ij1 == 3     .OR. ij2 == 3     .OR. ij3 == 3    ) )   lsend_bdyint(ib_bdy,igrd,3) = .true. 
     1209               IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. & 
     1210                    & ( ijbi == jpj-2 .OR. ij1 == jpj-2 .OR. ij2 == jpj-2 .OR. ij3 == jpj-2) )   lsend_bdyint(ib_bdy,igrd,4) = .true. 
    12131211               IF( ij == 2     .AND. ( nbondj ==  1 .OR. nbondj == 0 ) .AND. ijbe == 3     )   lsend_bdyext(ib_bdy,igrd,3) = .true. 
    12141212               IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. ijbe == jpj-2 )   lsend_bdyext(ib_bdy,igrd,4) = .true. 
     
    12341232         END IF 
    12351233      END DO 
    1236        
     1234 
    12371235      ! 
    12381236      ! Tidy up 
     
    12411239      ! 
    12421240   END SUBROUTINE bdy_segs 
     1241 
     1242 
     1243   SUBROUTINE find_neib( ii, ij, itreat, ii1, ij1, ii2, ij2, ii3, ij3 ) 
     1244     !!---------------------------------------------------------------------- 
     1245     !!                 ***  ROUTINE find_neib  *** 
     1246     !! 
     1247     !! ** Purpose :   get ii1, ij1, ii2, ij2, ii3, ij3, the indices of 
     1248     !!               the free ocean neighbours of (ii,ij) for bdy treatment 
     1249     !! 
     1250     !!---------------------------------------------------------------------- 
     1251     INTEGER, INTENT(in   )      ::   ii, ij, itreat 
     1252     INTEGER, INTENT(  out)      ::   ii1, ij1, ii2, ij2, ii3, ij3 
     1253     !!---------------------------------------------------------------------- 
     1254     SELECT CASE( itreat )   ! points that will be used by bdy routines, -1 will be discarded 
     1255        !               !               !     _____     !     _____      
     1256        !  1 |   o      !  2  o   |     !  3 | x        !  4     x |     
     1257        !    |_x_ _     !    _ _x_|     !    |   o      !      o   | 
     1258     CASE( 1 )    ;   ii1 = ii+1   ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1 
     1259     CASE( 2 )    ;   ii1 = ii-1   ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1 
     1260     CASE( 3 )    ;   ii1 = ii+1   ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1 
     1261     CASE( 4 )    ;   ii1 = ii-1   ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1 
     1262        !    |          !         |     !      o        !    ______                   ! or incomplete corner 
     1263        ! 5  | x o      ! 6   o x |     ! 7  __x__      ! 8    x                      !  7  ____ o 
     1264        !    |          !         |     !               !      o                      !         |x___ 
     1265     CASE( 5 )    ;   ii1 = ii+1   ;   ij1 = ij     ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1 
     1266     CASE( 6 )    ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1 
     1267     CASE( 7 )    ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1 
     1268     CASE( 8 )    ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1 
     1269        !        o      !        o      !    _____|     !       |_____   
     1270        !  9 ____x o    ! 10   o x___   ! 11     x o    ! 12   o x       
     1271        !         |     !       |       !        o      !        o       
     1272     CASE( 9  )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1  
     1273     CASE( 10 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1 
     1274     CASE( 11 )   ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1 
     1275     CASE( 12 )   ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1 
     1276        !    |_  o      !        o  _|  !     ¨¨|_|¨¨   !       o          
     1277        ! 13  _| x o    !  14  o x |_   !  15  o x o    ! 16  o x o        
     1278        !    |   o      !        o   |  !        o      !    __|¨|__  
     1279     CASE( 13 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = ii     ;   ij3 = ij-1    
     1280     CASE( 14 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = ii     ;   ij3 = ij-1  
     1281     CASE( 15 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij-1   ;   ii3 = ii+1   ;   ij3 = ij    
     1282     CASE( 16 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij+1   ;   ii3 = ii+1   ;   ij3 = ij 
     1283     END SELECT 
     1284   END SUBROUTINE find_neib 
    12431285 
    12441286 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdylib.F90

    r11071 r11183  
    1515   USE bdy_oce        ! ocean open boundary conditions 
    1616   USE phycst         ! physical constants 
     17   USE bdyini 
    1718   ! 
    1819   USE in_out_manager ! 
     
    452453      INTEGER  ::   ii, ij   ! 2D addresses 
    453454      INTEGER  ::   ipkm1    ! size of phia third dimension minus 1 
    454       INTEGER  ::   flagu, flagv                        ! short cuts 
    455       INTEGER  ::   ii1, ii2, ii3, ij1, ij2, ij3 
     455      INTEGER  ::   ii1, ii2, ii3, ij1, ij2, ij3, itreat 
    456456      !!---------------------------------------------------------------------- 
    457457      ! 
     
    468468         ii = idx%nbi(ib,igrd) 
    469469         ij = idx%nbj(ib,igrd) 
    470          ! 
    471          SELECT CASE( idx%ntreat(ib,igrd) )   ! select free ocean neighbours 
    472             !     o  
    473             !  ___x___         ! either flagu or flagv = 0 
    474             CASE( 0 )   ;   flagu = NINT( idx%flagu(ib,igrd) )   ;   ii1 = ii+flagu 
    475                             flagv = NINT( idx%flagv(ib,igrd) )   ;   ij1 = ij+flagv 
    476             !               !              !     _____     !     _____      
    477             !  1 |   o      !  2  o   |    !  3 | x        !  4     x |     
    478             !    |_x_ _     !    _ _x_|    !    |   o      !      o   | 
    479             CASE( 1 )   ;   ii1 = ii+1   ;   ij1 = ij+1 
    480             CASE( 2 )   ;   ii1 = ii-1   ;   ij1 = ij+1 
    481             CASE( 3 )   ;   ii1 = ii+1   ;   ij1 = ij-1 
    482             CASE( 4 )   ;   ii1 = ii-1   ;   ij1 = ij-1 
    483             !        o      !        o      !    _____|       !       |_____   
    484             !  5 ____x o    !  6   o x___   ! 7      x o      !  8   o x       
    485             !         |     !       |       !        o        !        o       
    486             CASE( 5 )   ;   ii1 = ii   ;   ij1 = ij+1   ;   ii2 = ii+1   ;   ij2 = ij    
    487             CASE( 6 )   ;   ii1 = ii   ;   ij1 = ij+1   ;   ii2 = ii-1   ;   ij2 = ij 
    488             CASE( 7 )   ;   ii1 = ii   ;   ij1 = ij-1   ;   ii2 = ii+1   ;   ij2 = ij   
    489             CASE( 8 )   ;   ii1 = ii   ;   ij1 = ij-1   ;   ii2 = ii-1   ;   ij2 = ij   
    490             !    |_  o      !        o  _|  !     ¨¨|_|¨¨   !       o          
    491             !  9  _| x o    !  10  o x |_   !  11  o x o    ! 12  o x o        
    492             !    |   o      !        o   |  !        o      !    __|¨|__  
    493             CASE( 9  )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = ii     ;   ij3 = ij-1    
    494             CASE( 10 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = ii     ;   ij3 = ij-1  
    495             CASE( 11 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij-1   ;   ii3 = ii+1   ;   ij3 = ij    
    496             CASE( 12 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij+1   ;   ii3 = ii+1   ;   ij3 = ij 
     470         itreat = idx%ntreat(ib,igrd) 
     471         CALL find_neib( ii, ij, itreat, ii1, ij1, ii2, ij2, ii3, ij3 ) 
     472         SELECT CASE( itreat ) 
     473         CASE( 1:8 ) 
     474            IF( ii1 < 1 .OR. ii1 > jpi .OR. ij1 < 1 .OR. ij1 > jpj )   CYCLE 
     475            DO ik = 1, ipkm1 
     476               IF( pmask(ii1,ij1,ik) /= 0. )   phia(ii,ij,ik) = phia(ii1,ij1,ik)   
     477            END DO 
     478         CASE( 9:12 ) 
     479            IF( ii1 < 1 .OR. ii1 > jpi .OR. ij1 < 1 .OR. ij1 > jpj )   CYCLE 
     480            IF( ii2 < 1 .OR. ii2 > jpi .OR. ij2 < 1 .OR. ij2 > jpj )   CYCLE 
     481            DO ik = 1, ipkm1 
     482               zweight = pmask(ii1,ij1,ik) + pmask(ii2,ij2,ik) 
     483               IF( zweight /= 0. )   phia(ii,ij,ik) = ( phia(ii1,ij1,ik) + phia(ii2,ij2,ik) ) / zweight 
     484            END DO 
     485         CASE( 13:16 ) 
     486            IF( ii1 < 1 .OR. ii1 > jpi .OR. ij1 < 1 .OR. ij1 > jpj )   CYCLE 
     487            IF( ii2 < 1 .OR. ii2 > jpi .OR. ij2 < 1 .OR. ij2 > jpj )   CYCLE 
     488            IF( ii3 < 1 .OR. ii3 > jpi .OR. ij3 < 1 .OR. ij3 > jpj )   CYCLE 
     489            DO ik = 1, ipkm1 
     490               zweight = pmask(ii1,ij1,ik) + pmask(ii2,ij2,ik) + pmask(ii3,ij3,ik) 
     491               IF( zweight /= 0. )   phia(ii,ij,ik) = ( phia(ii1,ij1,ik) + phia(ii2,ij2,ik) + phia(ii3,ij3,ik) ) / zweight 
     492            END DO 
    497493         END SELECT 
    498          ! 
    499          SELECT CASE( idx%ntreat(ib,igrd) ) 
    500             CASE( 0:4 ) 
    501                IF( ii1 < 1 .OR. ii1 > jpi .OR. ij1 < 1 .OR. ij1 > jpj )   CYCLE 
    502                DO ik = 1, ipkm1 
    503                   IF( pmask(ii1,ij1,ik) /= 0. )   phia(ii,ij,ik) = phia(ii1,ij1,ik)   
    504                END DO 
    505             CASE( 5:8 ) 
    506                IF( ii1 < 1 .OR. ii1 > jpi .OR. ij1 < 1 .OR. ij1 > jpj )   CYCLE 
    507                IF( ii2 < 1 .OR. ii2 > jpi .OR. ij2 < 1 .OR. ij2 > jpj )   CYCLE 
    508                DO ik = 1, ipkm1 
    509                   zweight = pmask(ii1,ij1,ik) + pmask(ii2,ij2,ik) 
    510                   IF( zweight /= 0. )   phia(ii,ij,ik) = ( phia(ii1,ij1,ik) + phia (ii2,ij2,ik) ) / zweight 
    511                END DO 
    512             CASE( 9:12 ) 
    513                IF( ii1 < 1 .OR. ii1 > jpi .OR. ij1 < 1 .OR. ij1 > jpj )   CYCLE 
    514                IF( ii2 < 1 .OR. ii2 > jpi .OR. ij2 < 1 .OR. ij2 > jpj )   CYCLE 
    515                IF( ii3 < 1 .OR. ii3 > jpi .OR. ij3 < 1 .OR. ij3 > jpj )   CYCLE 
    516                DO ik = 1, ipkm1 
    517                   zweight = pmask(ii1,ij1,ik) + pmask(ii2,ij2,ik) + pmask(ii3,ij3,ik) 
    518                   IF( zweight /= 0. )   phia(ii,ij,ik) = ( phia(ii1,ij1,ik) + phia (ii2,ij2,ik) + phia (ii3,ij3,ik) ) / zweight 
    519                END DO 
    520          END SELECT 
    521          ! 
    522494      END DO 
    523495      ! 
Note: See TracChangeset for help on using the changeset viewer.