Changeset 10679


Ignore:
Timestamp:
2019-02-14T14:11:43+01:00 (18 months ago)
Author:
mathiot
Message:

branch for solution 2 of ticket #2238

Location:
NEMO/branches/2019/fix_ticket2238_solution2
Files:
7 edited
1 copied

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB/icb_oce.F90

    r10425 r10679  
    6363      REAL(wp) ::   uo, vo, ui, vi, ua, va, ssh_x, ssh_y, sst, cn, hi    ! properties of iceberg environment  
    6464      REAL(wp) ::   mass_of_bits, heat_density 
     65      REAL(wp), DIMENSION(4) :: xRK1, xRK2, xRK3, xRK4 
     66      REAL(wp), DIMENSION(4) :: yRK1, yRK2, yRK3, yRK4 
    6567   END TYPE point 
    6668 
     
    8587   ! Extra arrays with bigger halo, needed when interpolating forcing onto iceberg position 
    8688   ! particularly for MPP when iceberg can lie inside T grid but outside U, V, or f grid 
    87    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   uo_e, vo_e 
    88    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ff_e, tt_e, fr_e, hicth 
    89    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ua_e, va_e 
    90    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ssh_e 
    91 #if defined key_si3 || defined key_cice 
    92    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ui_e, vi_e 
    93 #endif 
     89   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hicth 
    9490 
    9591   !!gm almost all those PARAM ARE defined in NEMO 
     
    168164      icb_alloc = icb_alloc + ill 
    169165      ! 
    170       ! expanded arrays for bilinear interpolation 
    171       ALLOCATE( uo_e(0:jpi+1,0:jpj+1) , ua_e(0:jpi+1,0:jpj+1) ,   & 
    172          &      vo_e(0:jpi+1,0:jpj+1) , va_e(0:jpi+1,0:jpj+1) ,   & 
    173 #if defined key_si3 || defined key_cice 
    174          &      ui_e(0:jpi+1,0:jpj+1) ,                            & 
    175          &      vi_e(0:jpi+1,0:jpj+1) ,                            & 
    176 #endif 
    177          &      ff_e(0:jpi+1,0:jpj+1) , fr_e(0:jpi+1,0:jpj+1)  ,   & 
    178          &      tt_e(0:jpi+1,0:jpj+1) , ssh_e(0:jpi+1,0:jpj+1) ,   & 
    179          &      hicth(0:jpi+1,0:jpj+1),                            & 
    180          &      first_width(nclasses) , first_length(nclasses) ,   & 
     166      ALLOCATE( first_width(nclasses) , first_length(nclasses) ,   & 
    181167         &      src_calving (jpi,jpj) ,                            & 
    182168         &      src_calving_hflx(jpi,jpj) , STAT=ill) 
  • NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB/icbdyn.F90

    r10570 r10679  
    1818   USE icbutl         ! iceberg utility routines 
    1919   USE icbdia         ! iceberg budget routines 
     20   USE icblbc         ! iceberg: lateral boundary routines (including mpp) 
     21   ! 
     22   USE in_out_manager 
     23   USE lib_mpp        ! massively parallel library  
    2024 
    2125   IMPLICIT NONE 
     
    6872      zdt_6 = zdt / 6._wp 
    6973 
    70       berg => first_berg                    ! start from the first berg 
    71       ! 
    72       DO WHILE ( ASSOCIATED(berg) )          !==  loop over all bergs  ==! 
     74      IF( ASSOCIATED(first_berg) ) THEN 
     75         berg => first_berg                    ! start from the first berg 
    7376         ! 
    74          pt => berg%current_point 
    75  
    76          ll_bounced = .FALSE. 
    77  
    78  
    79          ! STEP 1 ! 
    80          ! ====== ! 
    81          zxi1 = pt%xi   ;   zuvel1 = pt%uvel     !**   X1 in (i,j)  ;  V1 in m/s 
    82          zyj1 = pt%yj   ;   zvvel1 = pt%vvel 
    83  
    84  
    85          !                                         !**   A1 = A(X1,V1) 
    86          CALL icb_accel( berg , zxi1, ze1, zuvel1, zuvel1, zax1,     & 
    87             &                   zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2 ) 
    88          ! 
    89          zu1 = zuvel1 / ze1                           !**   V1 in d(i,j)/dt 
    90          zv1 = zvvel1 / ze2 
    91  
    92          ! STEP 2 ! 
    93          ! ====== ! 
    94          !                                         !**   X2 = X1+dt/2*V1   ;   V2 = V1+dt/2*A1 
    95          ! position using di/dt & djdt   !   V2  in m/s 
    96          zxi2 = zxi1 + zdt_2 * zu1          ;   zuvel2 = zuvel1 + zdt_2 * zax1 
    97          zyj2 = zyj1 + zdt_2 * zv1          ;   zvvel2 = zvvel1 + zdt_2 * zay1 
    98          ! 
    99          CALL icb_ground( zxi2, zxi1, zu1,   & 
    100             &             zyj2, zyj1, zv1, ll_bounced ) 
    101  
    102          !                                         !**   A2 = A(X2,V2) 
    103          CALL icb_accel( berg , zxi2, ze1, zuvel2, zuvel1, zax2,    & 
    104             &                   zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 ) 
    105          ! 
    106          zu2 = zuvel2 / ze1                           !**   V2 in d(i,j)/dt 
    107          zv2 = zvvel2 / ze2 
    108          ! 
    109          ! STEP 3 ! 
    110          ! ====== ! 
    111          !                                         !**  X3 = X1+dt/2*V2  ;   V3 = V1+dt/2*A2; A3=A(X3) 
    112          zxi3  = zxi1  + zdt_2 * zu2   ;   zuvel3 = zuvel1 + zdt_2 * zax2 
    113          zyj3  = zyj1  + zdt_2 * zv2   ;   zvvel3 = zvvel1 + zdt_2 * zay2 
    114          ! 
    115          CALL icb_ground( zxi3, zxi1, zu3,   & 
    116             &                zyj3, zyj1, zv3, ll_bounced ) 
    117  
    118          !                                         !**   A3 = A(X3,V3) 
    119          CALL icb_accel( berg , zxi3, ze1, zuvel3, zuvel1, zax3,    & 
    120             &                   zyj3, ze2, zvvel3, zvvel1, zay3, zdt ) 
    121          ! 
    122          zu3 = zuvel3 / ze1                           !**   V3 in d(i,j)/dt 
    123          zv3 = zvvel3 / ze2 
    124  
    125          ! STEP 4 ! 
    126          ! ====== ! 
    127          !                                         !**   X4 = X1+dt*V3   ;   V4 = V1+dt*A3 
    128          zxi4 = zxi1 + zdt * zu3   ;   zuvel4 = zuvel1 + zdt * zax3 
    129          zyj4 = zyj1 + zdt * zv3   ;   zvvel4 = zvvel1 + zdt * zay3 
    130  
    131          CALL icb_ground( zxi4, zxi1, zu4,   & 
    132             &             zyj4, zyj1, zv4, ll_bounced ) 
    133  
    134          !                                         !**   A4 = A(X4,V4) 
    135          CALL icb_accel( berg , zxi4, ze1, zuvel4, zuvel1, zax4,    & 
    136             &                   zyj4, ze2, zvvel4, zvvel1, zay4, zdt ) 
    137  
    138          zu4 = zuvel4 / ze1                           !**   V4 in d(i,j)/dt 
    139          zv4 = zvvel4 / ze2 
    140  
    141          ! FINAL STEP ! 
    142          ! ========== ! 
    143          !                                         !**   Xn = X1+dt*(V1+2*V2+2*V3+V4)/6 
    144          !                                         !**   Vn = V1+dt*(A1+2*A2+2*A3+A4)/6 
    145          zxi_n   = pt%xi   + zdt_6 * (  zu1  + 2.*(zu2  + zu3 ) + zu4  ) 
    146          zyj_n   = pt%yj   + zdt_6 * (  zv1  + 2.*(zv2  + zv3 ) + zv4  ) 
    147          zuvel_n = pt%uvel + zdt_6 * (  zax1 + 2.*(zax2 + zax3) + zax4 ) 
    148          zvvel_n = pt%vvel + zdt_6 * (  zay1 + 2.*(zay2 + zay3) + zay4 ) 
    149  
    150          CALL icb_ground( zxi_n, zxi1, zuvel_n,   & 
    151             &             zyj_n, zyj1, zvvel_n, ll_bounced ) 
    152  
    153          pt%uvel = zuvel_n                        !** save in berg structure 
    154          pt%vvel = zvvel_n 
    155          pt%xi   = zxi_n 
    156          pt%yj   = zyj_n 
    157  
    158          ! update actual position 
    159          pt%lon  = icb_utl_bilin_x(glamt, pt%xi, pt%yj ) 
    160          pt%lat  = icb_utl_bilin(gphit, pt%xi, pt%yj, 'T' ) 
    161  
    162          berg => berg%next                         ! switch to the next berg 
    163          ! 
    164       END DO                                  !==  end loop over all bergs  ==! 
     77         DO WHILE ( ASSOCIATED(berg) )          !==  loop over all bergs  ==! 
     78            ! 
     79            pt => berg%current_point 
     80            ! 
     81            ! STEP 1 ! 
     82            ! ====== ! 
     83            zxi1 = pt%xi   ;   zuvel1 = pt%uvel     !**   X1 in (i,j)  ;  V1 in m/s 
     84            zyj1 = pt%yj   ;   zvvel1 = pt%vvel 
     85            ! 
     86            !                                         !**   A1 = A(X1,V1) 
     87            CALL icb_accel( berg , zxi1, ze1, zuvel1, zuvel1, zax1,     & 
     88               &                   zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2 ) 
     89            ! 
     90            zu1 = zuvel1 / ze1                           !**   V1 in d(i,j)/dt 
     91            zv1 = zvvel1 / ze2 
     92            ! 
     93            ! STEP 2 ! 
     94            ! ====== ! 
     95            !                                         !**   X2 = X1+dt/2*V1   ;   V2 = V1+dt/2*A1 
     96            ! position using di/dt & djdt   !   V2  in m/s 
     97            zxi2 = zxi1 + zdt_2 * zu1          ;   zuvel2 = zuvel1 + zdt_2 * zax1 
     98            zyj2 = zyj1 + zdt_2 * zv1          ;   zvvel2 = zvvel1 + zdt_2 * zay1 
     99            ! 
     100            ! check ground iceberg 
     101            CALL icb_ground( zxi2, zxi1, zu1,   & 
     102               &             zyj2, zyj1, zv1, ll_bounced ) 
     103            ! 
     104            ! SAVE intermediaite properties 
     105            pt%xRK1 = (/ zxi1, zuvel1, zu1, zax1 /); pt%yRK1 = (/ zyj1, zvvel1, zv1, zay1 /);  
     106            pt%xRK2 = (/ zxi2, zuvel2, 0.0, 0.0  /); pt%yRK2 = (/ zyj2, zvvel2, 0.0, 0.0  /);  
     107            pt%xi = zxi2 ; pt%yj = zyj2 ; 
     108            !  
     109            berg => berg%next                         ! switch to the next berg 
     110            ! 
     111         END DO 
     112      END IF 
     113      ! 
     114      ! EXCHANGE BERG (in case first step fall out of the inner domain) 
     115      IF( lk_mpp ) THEN   ;          CALL icb_lbc_mpp()       ! Send bergs to other PEs 
     116      ELSE                ;          CALL icb_lbc()           ! Deal with any cyclic boundaries in non-mpp case 
     117      ENDIF 
     118      ! 
     119      IF( ASSOCIATED(first_berg) ) THEN 
     120         berg => first_berg                    ! start from the first berg 
     121         DO WHILE ( ASSOCIATED(berg) )          !==  loop over all bergs  ==! 
     122            ! 
     123            pt => berg%current_point 
     124            !  
     125            ! LOAD intermediaite properties 
     126            zxi1 = pt%xRK1(1) ; zuvel1 = pt%xRK1(2) ; zu1 = pt%xRK1(3) ; zax1 = pt%xRK1(4)  
     127            zxi2 = pt%xRK2(1) ; zuvel2 = pt%xRK2(2) ; zu2 = pt%xRK2(3) ; zax2 = pt%xRK2(4)  
     128            zyj1 = pt%yRK1(1) ; zvvel1 = pt%yRK1(2) ; zv1 = pt%yRK1(3) ; zay1 = pt%yRK1(4)  
     129            zyj2 = pt%yRK2(1) ; zvvel2 = pt%yRK2(2) ; zv2 = pt%yRK2(3) ; zay2 = pt%yRK2(4)  
     130            ! 
     131            !                                         !**   A2 = A(X2,V2) 
     132            CALL icb_accel( berg , zxi2, ze1, zuvel2, zuvel1, zax2,    & 
     133               &                   zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 ) 
     134            ! 
     135            zu2 = zuvel2 / ze1                           !**   V2 in d(i,j)/dt 
     136            zv2 = zvvel2 / ze2 
     137            ! 
     138            ! STEP 3 ! 
     139            ! ====== ! 
     140            !                                         !**  X3 = X1+dt/2*V2  ;   V3 = V1+dt/2*A2; A3=A(X3) 
     141            zxi3  = zxi1  + zdt_2 * zu2   ;   zuvel3 = zuvel1 + zdt_2 * zax2 
     142            zyj3  = zyj1  + zdt_2 * zv2   ;   zvvel3 = zvvel1 + zdt_2 * zay2 
     143            ! 
     144            CALL icb_ground( zxi3, zxi1, zu3,   & 
     145               &             zyj3, zyj1, zv3, ll_bounced ) 
     146            ! 
     147            ! SAVE intermediaite properties 
     148            pt%xRK2(3:4) =               (/ zu2, zax2 /); pt%yRK2(3:4) =               (/ zv2, zay2 /); 
     149            pt%xRK3      = (/ zxi3, zuvel3, 0.0, 0.0  /); pt%yRK3      = (/ zyj3, zvvel3, 0.0, 0.0  /);  
     150            pt%xi = zxi3 ; pt%yj = zyj3 ; 
     151            !  
     152            berg => berg%next                         ! switch to the next berg 
     153            ! 
     154         END DO 
     155      END IF 
     156      ! 
     157      ! EXCHANGE BERG (in case first step fall out of the inner domain) 
     158      IF( lk_mpp ) THEN   ;          CALL icb_lbc_mpp()       ! Send bergs to other PEs 
     159      ELSE                ;          CALL icb_lbc()           ! Deal with any cyclic boundaries in non-mpp case 
     160      ENDIF 
     161      ! 
     162      IF( ASSOCIATED(first_berg) ) THEN 
     163         berg => first_berg                    ! start from the first berg 
     164         DO WHILE ( ASSOCIATED(berg) )          !==  loop over all bergs  ==! 
     165            ! 
     166            pt => berg%current_point 
     167            !  
     168            ! LOAD intermediaite properties 
     169            zxi1 = pt%xRK1(1) ; zuvel1 = pt%xRK1(2) ; zu1 = pt%xRK1(3) ; zax1 = pt%xRK1(4)  
     170            zxi3 = pt%xRK3(1) ; zuvel3 = pt%xRK3(2) ; zu3 = pt%xRK3(3) ; zax3 = pt%xRK3(4)  
     171            zyj1 = pt%yRK1(1) ; zvvel1 = pt%yRK1(2) ; zv1 = pt%yRK1(3) ; zay1 = pt%yRK1(4)  
     172            zyj3 = pt%yRK3(1) ; zvvel3 = pt%yRK3(2) ; zv3 = pt%yRK3(3) ; zay3 = pt%yRK3(4)  
     173            ! 
     174            !                                         !**   A3 = A(X3,V3) 
     175            CALL icb_accel( berg , zxi3, ze1, zuvel3, zuvel1, zax3,    & 
     176               &                   zyj3, ze2, zvvel3, zvvel1, zay3, zdt ) 
     177            ! 
     178            zu3 = zuvel3 / ze1                        !**   V3 in d(i,j)/dt 
     179            zv3 = zvvel3 / ze2 
     180            ! 
     181            ! STEP 4 ! 
     182            ! ====== ! 
     183            !                                         !**   X4 = X1+dt*V3   ;   V4 = V1+dt*A3 
     184            zxi4 = zxi1 + zdt * zu3   ;   zuvel4 = zuvel1 + zdt * zax3 
     185            zyj4 = zyj1 + zdt * zv3   ;   zvvel4 = zvvel1 + zdt * zay3 
     186            ! 
     187            CALL icb_ground( zxi4, zxi1, zu4,   & 
     188               &             zyj4, zyj1, zv4, ll_bounced ) 
     189            ! 
     190            ! SAVE intermediaite properties 
     191            pt%xRK3(3:4) =               (/ zu3, zax3 /); pt%yRK3(3:4) =               (/ zv3, zay3 /);  
     192            pt%xRK4      = (/ zxi4, zuvel4, 0.0, 0.0  /); pt%yRK4      = (/ zyj4, zvvel4, 0.0, 0.0  /);  
     193            pt%xi = zxi4 ; pt%yj = zyj4 ; 
     194            ! 
     195            berg => berg%next                         ! switch to the next berg 
     196            ! 
     197         END DO 
     198      END IF 
     199      ! 
     200      ! EXCHANGE BERG (in case first step fall out of the inner domain) 
     201      IF( lk_mpp ) THEN   ;          CALL icb_lbc_mpp()       ! Send bergs to other PEs 
     202      ELSE                ;          CALL icb_lbc()           ! Deal with any cyclic boundaries in non-mpp case 
     203      ENDIF 
     204      ! 
     205      IF( ASSOCIATED(first_berg) ) THEN 
     206         berg => first_berg                    ! start from the first berg 
     207         DO WHILE ( ASSOCIATED(berg) )          !==  loop over all bergs  ==! 
     208            ! 
     209            pt => berg%current_point 
     210            !  
     211            ! LOAD intermediaite properties 
     212            zxi1 = pt%xRK1(1) ; zuvel1 = pt%xRK1(2) ; zu1 = pt%xRK1(3) ; zax1 = pt%xRK1(4)  
     213            zxi2 = pt%xRK2(1) ; zuvel2 = pt%xRK2(2) ; zu2 = pt%xRK2(3) ; zax2 = pt%xRK2(4)  
     214            zxi3 = pt%xRK3(1) ; zuvel3 = pt%xRK3(2) ; zu3 = pt%xRK3(3) ; zax3 = pt%xRK3(4)  
     215            zxi4 = pt%xRK4(1) ; zuvel4 = pt%xRK4(2) ; zu4 = pt%xRK4(3) ; zax4 = pt%xRK4(4)  
     216            zyj1 = pt%yRK1(1) ; zvvel1 = pt%yRK1(2) ; zv1 = pt%yRK1(3) ; zay1 = pt%yRK1(4)  
     217            zyj2 = pt%yRK2(1) ; zvvel2 = pt%yRK2(2) ; zv2 = pt%yRK2(3) ; zay2 = pt%yRK2(4)  
     218            zyj3 = pt%yRK3(1) ; zvvel3 = pt%yRK3(2) ; zv3 = pt%yRK3(3) ; zay3 = pt%yRK3(4)  
     219            zyj4 = pt%yRK4(1) ; zvvel4 = pt%yRK4(2) ; zv4 = pt%yRK4(3) ; zay4 = pt%yRK4(4)  
     220            ! 
     221            !                                         !**   A4 = A(X4,V4) 
     222            CALL icb_accel( berg , zxi4, ze1, zuvel4, zuvel1, zax4,    & 
     223               &                   zyj4, ze2, zvvel4, zvvel1, zay4, zdt ) 
     224            ! 
     225            zu4 = zuvel4 / ze1                        !**   V4 in d(i,j)/dt 
     226            zv4 = zvvel4 / ze2 
     227            ! 
     228            ! FINAL STEP ! 
     229            ! ========== ! 
     230            !                                         !**   Xn = X1+dt*(V1+2*V2+2*V3+V4)/6 
     231            !                                         !**   Vn = V1+dt*(A1+2*A2+2*A3+A4)/6 
     232            zxi_n   = zxi1   + zdt_6 * (  zu1  + 2._wp*(zu2  + zu3 ) + zu4  ) 
     233            zyj_n   = zyj1   + zdt_6 * (  zv1  + 2._wp*(zv2  + zv3 ) + zv4  ) 
     234            zuvel_n = zuvel1 + zdt_6 * (  zax1 + 2._wp*(zax2 + zax3) + zax4 ) 
     235            zvvel_n = zvvel1 + zdt_6 * (  zay1 + 2._wp*(zay2 + zay3) + zay4 ) 
     236            ! 
     237            CALL icb_ground( zxi_n, zxi1, zuvel_n,   & 
     238               &             zyj_n, zyj1, zvvel_n, ll_bounced ) 
     239            ! 
     240            !                                        !** save in berg structure 
     241            pt%xi   = zxi_n   ; pt%yj   = zyj_n ; 
     242            pt%uvel = zuvel_n ; pt%vvel = zvvel_n 
     243            ! 
     244            berg => berg%next                        ! switch to the next berg 
     245            ! 
     246         END DO 
     247      END IF 
     248      ! 
     249      ! EXCHANGE BERG (in case first step fall out of the inner domain) 
     250      IF( lk_mpp ) THEN   ;          CALL icb_lbc_mpp()       ! Send bergs to other PEs 
     251      ELSE                ;          CALL icb_lbc()           ! Deal with any cyclic boundaries in non-mpp case 
     252      ENDIF 
     253      ! 
     254      IF( ASSOCIATED(first_berg) ) THEN 
     255         berg => first_berg                    ! start from the first berg 
     256         DO WHILE ( ASSOCIATED(berg) )          !==  loop over all bergs  ==! 
     257            ! 
     258            pt => berg%current_point 
     259            ! 
     260            ! update actual position 
     261            pt%lon  = icb_utl_bilin_x(glamt, pt%xi, pt%yj ) 
     262            pt%lat  = icb_utl_bilin(gphit, pt%xi, pt%yj, 'T' ) 
     263            ! 
     264            berg => berg%next                         ! switch to the next berg 
     265            ! 
     266         END DO                                  !==  end loop over all bergs  ==! 
     267      END IF 
    165268      ! 
    166269   END SUBROUTINE icb_dyn 
     
    182285      LOGICAL , INTENT(  out) ::   ld_bounced   ! bounced indicator 
    183286      ! 
    184       INTEGER  ::   ii, ii0 
    185       INTEGER  ::   ij, ij0 
     287      INTEGER  ::   ii, iig, iig0 
     288      INTEGER  ::   ij, ijg, ijg0 
    186289      INTEGER  ::   ibounce_method 
    187290      !!---------------------------------------------------------------------- 
    188291      ! 
     292      ! (PM) is this variable used somewhere ? 
    189293      ld_bounced = .FALSE. 
    190294      ! 
    191       ii0 = INT( pi0+0.5 )   ;   ij0 = INT( pj0+0.5 )       ! initial gridpoint position (T-cell) 
    192       ii  = INT( pi +0.5 )   ;   ij  = INT( pj +0.5 )       ! current     -         - 
    193       ! 
    194       IF( ii == ii0  .AND.  ij == ij0  )   RETURN           ! berg remains in the same cell 
     295      ! work with global position in case pi and pi0 are on different domains 
     296      iig0 = INT( pi0 + 0.5_wp )   ;   ijg0 = INT( pj0 + 0.5_wp )       ! initial gridpoint position (T-cell) 
     297      iig  = INT( pi  + 0.5_wp )   ;   ijg  = INT( pj  + 0.5_wp )       ! current     -         - 
     298      ! 
     299      IF( iig == iig0  .AND.  ijg == ijg0  )   RETURN           ! berg remains in the same cell 
    195300      ! 
    196301      ! map into current processor 
    197       ii0 = mi1( ii0 ) 
    198       ij0 = mj1( ij0 ) 
    199       ii  = mi1( ii  ) 
    200       ij  = mj1( ij  ) 
     302      ii  = mi1( iig  ) 
     303      ij  = mj1( ijg  ) 
    201304      ! 
    202305      IF(  tmask(ii,ij,1)  /=   0._wp  )   RETURN           ! berg reach a new t-cell, but an ocean one 
     
    221324         pv = 0._wp 
    222325      CASE ( 2 ) 
    223          IF( ii0 /= ii ) THEN 
     326         IF( iig0 /= iig ) THEN 
    224327            pi = pi0                   ! return back to the initial position 
    225328            pu = 0._wp                 ! zeroing of velocity in the direction of the grounding 
    226329         ENDIF 
    227          IF( ij0 /= ij ) THEN 
     330         IF( ijg0 /= ijg ) THEN 
    228331            pj = pj0                   ! return back to the initial position 
    229332            pv = 0._wp                 ! zeroing of velocity in the direction of the grounding 
  • NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB/icblbc.F90

    r10570 r10679  
    5656 
    5757   INTEGER, PARAMETER, PRIVATE ::   jp_delta_buf = 25             ! Size by which to increment buffers 
    58    INTEGER, PARAMETER, PRIVATE ::   jp_buffer_width = 15+nkounts  ! items to store for each berg 
     58   INTEGER, PARAMETER, PRIVATE ::   jp_buffer_width = 47+nkounts  ! items to store for each berg 
    5959 
    6060#endif 
     
    9090         DO WHILE( ASSOCIATED(this) ) 
    9191            pt => this%current_point 
    92             iine = INT( pt%xi + 0.5 ) 
     92            iine = INT( pt%xi + 0.5_wp ) 
    9393            IF( iine > mig(nicbei) ) THEN 
    9494               pt%xi = ricb_right + MOD(pt%xi, 1._wp ) - 1._wp 
     
    125125      DO WHILE( ASSOCIATED(this) ) 
    126126         pt => this%current_point 
    127          ijne = INT( pt%yj + 0.5 ) 
     127         ijne = INT( pt%yj + 0.5_wp ) 
    128128         IF( ijne .GT. mjg(nicbej) ) THEN 
    129129            ! 
    130             iine = INT( pt%xi + 0.5 ) 
     130            iine = INT( pt%xi + 0.5_wp ) 
    131131            ipts  = nicbfldpts (mi1(iine)) 
    132132            ! 
     
    232232         DO WHILE (ASSOCIATED(this)) 
    233233            pt => this%current_point 
    234             iine = INT( pt%xi + 0.5 ) 
     234            iine = INT( pt%xi + 0.5_wp ) 
    235235            IF( ipe_E >= 0 .AND. iine > mig(nicbei) ) THEN 
    236236               tmpberg => this 
     
    242242               ENDIF 
    243243               ! deal with periodic case 
    244                tmpberg%current_point%xi = ricb_right + MOD(tmpberg%current_point%xi, 1._wp ) - 1._wp 
     244               tmpberg%current_point%xi      = ricb_right + MOD(tmpberg%current_point%xi, 1._wp ) - 1._wp 
     245               tmpberg%current_point%xRK1(1) = ricb_right + MOD(tmpberg%current_point%xRK1(1), 1._wp ) - 1._wp 
     246               tmpberg%current_point%xRK2(1) = ricb_right + MOD(tmpberg%current_point%xRK2(1), 1._wp ) - 1._wp 
     247               tmpberg%current_point%xRK3(1) = ricb_right + MOD(tmpberg%current_point%xRK3(1), 1._wp ) - 1._wp 
     248               tmpberg%current_point%xRK4(1) = ricb_right + MOD(tmpberg%current_point%xRK4(1), 1._wp ) - 1._wp 
    245249               ! now pack it into buffer and delete from list 
    246250               CALL icb_pack_into_buffer( tmpberg, obuffer_e, ibergs_to_send_e) 
     
    255259               ENDIF 
    256260               ! deal with periodic case 
    257                tmpberg%current_point%xi = ricb_left + MOD(tmpberg%current_point%xi, 1._wp ) 
     261               tmpberg%current_point%xi      = ricb_left + MOD(tmpberg%current_point%xi, 1._wp ) 
     262               tmpberg%current_point%xRK1(1) = ricb_left + MOD(tmpberg%current_point%xRK1(1), 1._wp ) 
     263               tmpberg%current_point%xRK2(1) = ricb_left + MOD(tmpberg%current_point%xRK2(1), 1._wp ) 
     264               tmpberg%current_point%xRK3(1) = ricb_left + MOD(tmpberg%current_point%xRK3(1), 1._wp ) 
     265               tmpberg%current_point%xRK4(1) = ricb_left + MOD(tmpberg%current_point%xRK4(1), 1._wp ) 
    258266               ! now pack it into buffer and delete from list 
    259267               CALL icb_pack_into_buffer( tmpberg, obuffer_w, ibergs_to_send_w) 
     
    370378         DO WHILE (ASSOCIATED(this)) 
    371379            pt => this%current_point 
    372             ijne = INT( pt%yj + 0.5 ) 
     380            ijne = INT( pt%yj + 0.5_wp ) 
    373381            IF( ipe_N >= 0 .AND. ijne .GT. mjg(nicbej) ) THEN 
    374382               tmpberg => this 
     
    537545         DO WHILE (ASSOCIATED(this)) 
    538546            pt => this%current_point 
    539             iine = INT( pt%xi + 0.5 ) 
    540             ijne = INT( pt%yj + 0.5 ) 
     547            iine = INT( pt%xi + 0.5_wp ) 
     548            ijne = INT( pt%yj + 0.5_wp ) 
    541549            IF( iine .LT. mig(nicbdi) .OR. & 
    542550                iine .GT. mig(nicbei) .OR. & 
     
    544552                ijne .GT. mjg(nicbej)) THEN 
    545553               i = i + 1 
    546                WRITE(numicb,*) 'berg lost in halo: ', this%number(:),iine,ijne 
    547                WRITE(numicb,*) '                   ', nimpp, njmpp 
    548                WRITE(numicb,*) '                   ', nicbdi, nicbei, nicbdj, nicbej 
     554               WRITE(numicb,*) 'berg lost in halo:  ', this%number(:),iine,ijne, pt%xi, pt%yj 
     555               WRITE(numicb,*) 'nimpp, njmpp        ', nimpp, njmpp 
     556               WRITE(numicb,*) 'nicb di, ei, dj, ej ', nicbdi     , nicbei     , nicbdj     , nicbej 
     557               WRITE(numicb,*) 'mijg di, ei, dj, ej ', mig(nicbdi), mig(nicbei), mjg(nicbdj), mjg(nicbej) 
    549558               CALL flush( numicb ) 
    550559            ENDIF 
     
    611620               DO WHILE (ASSOCIATED(this)) 
    612621                  pt => this%current_point 
    613                   iine = INT( pt%xi + 0.5 ) 
    614                   ijne = INT( pt%yj + 0.5 ) 
     622                  iine = INT( pt%xi + 0.5_wp ) 
     623                  ijne = INT( pt%yj + 0.5_wp ) 
    615624                  iproc = nicbflddest(mi1(iine)) 
    616625                  IF( ijne .GT. mjg(nicbej) ) THEN 
     
    691700               DO WHILE (ASSOCIATED(this)) 
    692701                  pt => this%current_point 
    693                   iine = INT( pt%xi + 0.5 ) 
    694                   ijne = INT( pt%yj + 0.5 ) 
     702                  iine = INT( pt%xi + 0.5_wp ) 
     703                  ijne = INT( pt%yj + 0.5_wp ) 
    695704                  ipts  = nicbfldpts (mi1(iine)) 
    696705                  iproc = nicbflddest(mi1(iine)) 
     
    791800      IF( .NOT. ASSOCIATED(pbuff) ) CALL icb_increase_buffer( pbuff, jp_delta_buf ) 
    792801      IF( kb .GT. pbuff%size ) CALL icb_increase_buffer( pbuff, jp_delta_buf ) 
     802      IF( kb .GT. pbuff%size ) PRINT *, 'SHOULD NOT SEE THIS' 
    793803 
    794804      !! pack points into buffer 
     
    809819      pbuff%data(14,kb) = berg%current_point%heat_density 
    810820 
    811       pbuff%data(15,kb) = berg%mass_scaling 
     821      pbuff%data(15:18,kb) = berg%current_point%xRK1 
     822      pbuff%data(19:22,kb) = berg%current_point%xRK2 
     823      pbuff%data(23:26,kb) = berg%current_point%xRK3 
     824      pbuff%data(27:30,kb) = berg%current_point%xRK4 
     825 
     826      pbuff%data(31:34,kb) = berg%current_point%yRK1 
     827      pbuff%data(35:38,kb) = berg%current_point%yRK2 
     828      pbuff%data(39:42,kb) = berg%current_point%yRK3 
     829      pbuff%data(43:46,kb) = berg%current_point%yRK4 
     830 
     831      pbuff%data(47,kb) = berg%mass_scaling 
    812832      DO k=1,nkounts 
    813          pbuff%data(15+k,kb) = REAL( berg%number(k), wp ) 
     833         pbuff%data(47+k,kb) = REAL( berg%number(k), wp ) 
    814834      END DO 
    815835      ! 
     
    844864      pt%heat_density   =      pbuff%data(14,kb) 
    845865 
    846       currentberg%mass_scaling =      pbuff%data(15,kb) 
     866      pt%xRK1 = pbuff%data(15:18,kb) 
     867      pt%xRK2 = pbuff%data(19:22,kb)  
     868      pt%xRK3 = pbuff%data(23:26,kb) 
     869      pt%xRK4 = pbuff%data(27:30,kb) 
     870 
     871      pt%yRK1 = pbuff%data(31:34,kb) 
     872      pt%yRK2 = pbuff%data(35:38,kb) 
     873      pt%yRK3 = pbuff%data(39:42,kb) 
     874      pt%yRK4 = pbuff%data(43:46,kb) 
     875 
     876      currentberg%mass_scaling =      pbuff%data(47,kb) 
    847877      DO ik = 1, nkounts 
    848          currentberg%number(ik) = INT( pbuff%data(15+ik,kb) ) 
     878         currentberg%number(ik) = INT( pbuff%data(47+ik,kb) ) 
    849879      END DO 
    850880      ! 
  • NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB/icbstp.F90

    r10570 r10679  
    2929   USE icbclv         ! iceberg: calving routines 
    3030   USE icbthm         ! iceberg: thermodynamics routines 
    31    USE icblbc         ! iceberg: lateral boundary routines (including mpp) 
    3231   USE icbtrj         ! iceberg: trajectory I/O routines 
    3332   USE icbdia         ! iceberg: budget 
     
    9190 9100 FORMAT('kt= ',i8, ' day= ',i8,' secs=',i8) 
    9291      ! 
    93       !                                   !* copy nemo forcing arrays into iceberg versions with extra halo 
    94       CALL icb_utl_copy()                 ! only necessary for variables not on T points 
    95       ! 
    96       ! 
    9792      !                       !==  process icebergs  ==! 
    9893      !                              ! 
     
    10499      !                       !==  For each berg, evolve  ==! 
    105100      ! 
    106       IF( ASSOCIATED(first_berg) )   CALL icb_dyn( kt )       ! ice berg dynamics 
    107  
    108       IF( lk_mpp ) THEN   ;          CALL icb_lbc_mpp()       ! Send bergs to other PEs 
    109       ELSE                ;          CALL icb_lbc()           ! Deal with any cyclic boundaries in non-mpp case 
    110       ENDIF 
     101                                     CALL icb_dyn( kt )       ! ice berg dynamics 
    111102 
    112103      IF( ASSOCIATED(first_berg) )   CALL icb_thm( kt )       ! Ice berg thermodynamics (melting) + rolling 
  • NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB/icbutl.F90

    r10570 r10679  
    3131   PRIVATE 
    3232 
    33    PUBLIC   icb_utl_copy          ! routine called in icbstp module 
    3433   PUBLIC   icb_utl_interp        ! routine called in icbdyn, icbthm modules 
    3534   PUBLIC   icb_utl_bilin         ! routine called in icbini, icbdyn modules 
     
    5453CONTAINS 
    5554 
    56    SUBROUTINE icb_utl_copy() 
    57       !!---------------------------------------------------------------------- 
    58       !!                  ***  ROUTINE icb_utl_copy  *** 
    59       !! 
    60       !! ** Purpose :   iceberg initialization. 
    61       !! 
    62       !! ** Method  : - blah blah 
    63       !!---------------------------------------------------------------------- 
    64 #if defined key_si3 
    65       REAL(wp), DIMENSION(jpi,jpj) :: zssh_lead_m    !    ocean surface (ssh_m) if ice is not embedded 
    66       !                                              !    ocean surface in leads if ice is embedded    
    67 #endif 
    68       ! copy nemo forcing arrays into iceberg versions with extra halo 
    69       ! only necessary for variables not on T points 
    70       ! and ssh which is used to calculate gradients 
    71  
    72       uo_e(:,:) = 0._wp   ;   uo_e(1:jpi,1:jpj) = ssu_m(:,:) * umask(:,:,1) 
    73       vo_e(:,:) = 0._wp   ;   vo_e(1:jpi,1:jpj) = ssv_m(:,:) * vmask(:,:,1) 
    74       ff_e(:,:) = 0._wp   ;   ff_e(1:jpi,1:jpj) = ff_f (:,:)  
    75       tt_e(:,:) = 0._wp   ;   tt_e(1:jpi,1:jpj) = sst_m(:,:) 
    76       fr_e(:,:) = 0._wp   ;   fr_e(1:jpi,1:jpj) = fr_i (:,:) 
    77       ua_e(:,:) = 0._wp   ;   ua_e(1:jpi,1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk 
    78       va_e(:,:) = 0._wp   ;   va_e(1:jpi,1:jpj) = vtau (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk 
    79       ! 
    80       CALL lbc_lnk_icb( 'icbutl', uo_e, 'U', -1._wp, 1, 1 ) 
    81       CALL lbc_lnk_icb( 'icbutl', vo_e, 'V', -1._wp, 1, 1 ) 
    82       CALL lbc_lnk_icb( 'icbutl', ff_e, 'F', +1._wp, 1, 1 ) 
    83       CALL lbc_lnk_icb( 'icbutl', ua_e, 'U', -1._wp, 1, 1 ) 
    84       CALL lbc_lnk_icb( 'icbutl', va_e, 'V', -1._wp, 1, 1 ) 
    85       CALL lbc_lnk_icb( 'icbutl', fr_e, 'T', +1._wp, 1, 1 ) 
    86       CALL lbc_lnk_icb( 'icbutl', tt_e, 'T', +1._wp, 1, 1 ) 
    87 #if defined key_si3 
    88       hicth(:,:) = 0._wp ;  hicth(1:jpi,1:jpj) = hm_i (:,:)   
    89       ui_e(:,:) = 0._wp ;   ui_e(1:jpi, 1:jpj) = u_ice(:,:) 
    90       vi_e(:,:) = 0._wp ;   vi_e(1:jpi, 1:jpj) = v_ice(:,:) 
    91       !       
    92       ! compute ssh slope using ssh_lead if embedded 
    93       zssh_lead_m(:,:) = ice_var_sshdyn(ssh_m, snwice_mass, snwice_mass_b) 
    94       ssh_e(:,:) = 0._wp ;  ssh_e(1:jpi, 1:jpj) = zssh_lead_m(:,:) * tmask(:,:,1) 
    95       ! 
    96       CALL lbc_lnk_icb( 'icbutl', hicth, 'T', +1._wp, 1, 1 ) 
    97       CALL lbc_lnk_icb( 'icbutl', ui_e , 'U', -1._wp, 1, 1 ) 
    98       CALL lbc_lnk_icb( 'icbutl', vi_e , 'V', -1._wp, 1, 1 ) 
    99 #else 
    100       ssh_e(:,:) = 0._wp ;  ssh_e(1:jpi, 1:jpj) = ssh_m(:,:) * tmask(:,:,1) 
    101 #endif 
    102  
    103       !! special for ssh which is used to calculate slope 
    104       !! so fudge some numbers all the way around the boundary 
    105       ssh_e(0    ,    :) = ssh_e(1  ,  :) 
    106       ssh_e(jpi+1,    :) = ssh_e(jpi,  :) 
    107       ssh_e(:    ,    0) = ssh_e(:  ,  1) 
    108       ssh_e(:    ,jpj+1) = ssh_e(:  ,jpj) 
    109       ssh_e(0,0)         = ssh_e(1,1) 
    110       ssh_e(jpi+1,0)     = ssh_e(jpi,1) 
    111       ssh_e(0,jpj+1)     = ssh_e(1,jpj) 
    112       ssh_e(jpi+1,jpj+1) = ssh_e(jpi,jpj) 
    113       CALL lbc_lnk_icb( 'icbutl', ssh_e, 'T', +1._wp, 1, 1 ) 
    114       ! 
    115    END SUBROUTINE icb_utl_copy 
    116  
    117  
    11855   SUBROUTINE icb_utl_interp( pi, pe1, puo, pui, pua, pssh_i,   & 
    11956      &                       pj, pe2, pvo, pvi, pva, pssh_j,   & 
     
    14077      ! 
    14178      REAL(wp) ::   zcd, zmod       ! local scalars 
     79      REAL(wp), DIMENSION(jpi,jpj) :: zssh 
    14280      !!---------------------------------------------------------------------- 
    14381 
     
    14583      pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj ) 
    14684      ! 
    147       puo  = icb_utl_bilin_h( uo_e, pi, pj, 'U' )             ! ocean velocities 
    148       pvo  = icb_utl_bilin_h( vo_e, pi, pj, 'V' ) 
    149       psst = icb_utl_bilin_h( tt_e, pi, pj, 'T' )             ! SST 
    150       pcn  = icb_utl_bilin_h( fr_e , pi, pj, 'T' )            ! ice concentration 
    151       pff  = icb_utl_bilin_h( ff_e , pi, pj, 'F' )            ! Coriolis parameter 
    152       ! 
    153       pua  = icb_utl_bilin_h( ua_e , pi, pj, 'U' )            ! 10m wind 
    154       pva  = icb_utl_bilin_h( va_e , pi, pj, 'V' )            ! here (ua,va) are stress => rough conversion from stress to speed 
     85      puo  = icb_utl_bilin( ssu_m(:,:) * umask(:,:,1), pi, pj, 'U' )             ! ocean velocities 
     86      pvo  = icb_utl_bilin( ssv_m(:,:) * vmask(:,:,1), pi, pj, 'V' ) 
     87      psst = icb_utl_bilin( sst_m(:,:), pi, pj, 'T' )             ! SST 
     88      pcn  = icb_utl_bilin( fr_i(:,:) , pi, pj, 'T' )            ! ice concentration 
     89      pff  = icb_utl_bilin( ff_f(:,:) , pi, pj, 'F' )            ! Coriolis parameter 
     90      ! 
     91      pua  = icb_utl_bilin( utau(:,:) * umask(:,:,1), pi, pj, 'U' )            ! 10m wind 
     92      pva  = icb_utl_bilin( vtau(:,:) * vmask(:,:,1), pi, pj, 'V' )            ! here (ua,va) are stress => rough conversion from stress to speed 
    15593      zcd  = 1.22_wp * 1.5e-3_wp                              ! air density * drag coefficient 
    15694      zmod = 1._wp / MAX(  1.e-20, SQRT(  zcd * SQRT( pua*pua + pva*pva)  )  ) 
     
    15997 
    16098#if defined key_si3 
    161       pui = icb_utl_bilin_h( ui_e , pi, pj, 'U' )              ! sea-ice velocities 
    162       pvi = icb_utl_bilin_h( vi_e , pi, pj, 'V' ) 
    163       phi = icb_utl_bilin_h( hicth, pi, pj, 'T' )              ! ice thickness 
     99      pui = icb_utl_bilin( u_ice , pi, pj, 'U' )              ! sea-ice velocities 
     100      pvi = icb_utl_bilin( v_ice , pi, pj, 'V' ) 
     101      phi = icb_utl_bilin( hm_i  , pi, pj, 'T' )              ! ice thickness 
     102      ! Estimate SSH gradient in i- and j-direction (centred evaluation) 
     103      ! compute ssh slope using ssh_lead if embedded 
     104      zssh(:,:) = ice_var_sshdyn(ssh_m, snwice_mass, snwice_mass_b) 
    164105#else 
    165106      pui = 0._wp 
    166107      pvi = 0._wp 
    167108      phi = 0._wp 
     109      zssh(:,:) = ssh_m(:,:) 
    168110#endif 
     111      zssh(:,:) = zssh(:,:) * tmask(:,:,1) 
    169112 
    170113      ! Estimate SSH gradient in i- and j-direction (centred evaluation) 
    171       pssh_i = ( icb_utl_bilin_h( ssh_e, pi+0.1_wp, pj, 'T' ) -   & 
    172          &       icb_utl_bilin_h( ssh_e, pi-0.1_wp, pj, 'T' )  ) / ( 0.2_wp * pe1 ) 
    173       pssh_j = ( icb_utl_bilin_h( ssh_e, pi, pj+0.1_wp, 'T' ) -   & 
    174          &       icb_utl_bilin_h( ssh_e, pi, pj-0.1_wp, 'T' )  ) / ( 0.2_wp * pe2 ) 
     114      pssh_i = ( icb_utl_bilin( zssh(:,:), pi+0.1_wp, pj, 'T' ) -   & 
     115         &       icb_utl_bilin( zssh(:,:), pi-0.1_wp, pj, 'T' )  ) / ( 0.2_wp * pe1 ) 
     116      pssh_j = ( icb_utl_bilin( zssh(:,:), pi, pj+0.1_wp, 'T' ) -   & 
     117         &       icb_utl_bilin( zssh(:,:), pi, pj-0.1_wp, 'T' )  ) / ( 0.2_wp * pe2 ) 
    175118      ! 
    176119   END SUBROUTINE icb_utl_interp 
    177  
    178  
    179    REAL(wp) FUNCTION icb_utl_bilin_h( pfld, pi, pj, cd_type ) 
    180       !!---------------------------------------------------------------------- 
    181       !!                  ***  FUNCTION icb_utl_bilin  *** 
    182       !! 
    183       !! ** Purpose :   bilinear interpolation at berg location depending on the grid-point type 
    184       !!                this version deals with extra halo points 
    185       !! 
    186       !!       !!gm  CAUTION an optional argument should be added to handle 
    187       !!             the slip/no-slip conditions  ==>>> to be done later 
    188       !! 
    189       !!---------------------------------------------------------------------- 
    190       REAL(wp), DIMENSION(0:jpi+1,0:jpj+1), INTENT(in) ::   pfld      ! field to be interpolated 
    191       REAL(wp)                            , INTENT(in) ::   pi, pj    ! targeted coordinates in (i,j) referential 
    192       CHARACTER(len=1)                    , INTENT(in) ::   cd_type   ! type of pfld array grid-points: = T , U , V or F points 
    193       ! 
    194       INTEGER  ::   ii, ij   ! local integer 
    195       REAL(wp) ::   zi, zj   ! local real 
    196       !!---------------------------------------------------------------------- 
    197       ! 
    198       SELECT CASE ( cd_type ) 
    199       CASE ( 'T' ) 
    200          ! note that here there is no +0.5 added 
    201          ! since we're looking for four T points containing quadrant we're in of  
    202          ! current T cell 
    203          ii = MAX(1, INT( pi     )) 
    204          ij = MAX(1, INT( pj     ))    ! T-point 
    205          zi = pi - REAL(ii,wp) 
    206          zj = pj - REAL(ij,wp) 
    207       CASE ( 'U' ) 
    208          ii = MAX(1, INT( pi-0.5 )) 
    209          ij = MAX(1, INT( pj     ))    ! U-point 
    210          zi = pi - 0.5 - REAL(ii,wp) 
    211          zj = pj - REAL(ij,wp) 
    212       CASE ( 'V' ) 
    213          ii = MAX(1, INT( pi     )) 
    214          ij = MAX(1, INT( pj-0.5 ))    ! V-point 
    215          zi = pi - REAL(ii,wp) 
    216          zj = pj - 0.5 - REAL(ij,wp) 
    217       CASE ( 'F' ) 
    218          ii = MAX(1, INT( pi-0.5 )) 
    219          ij = MAX(1, INT( pj-0.5 ))    ! F-point 
    220          zi = pi - 0.5 - REAL(ii,wp) 
    221          zj = pj - 0.5 - REAL(ij,wp) 
    222       END SELECT 
    223       ! 
    224       ! find position in this processor. Prevent near edge problems (see #1389) 
    225       ! 
    226       IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1 
    227       ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi 
    228       ELSE                           ;   ii = mi1(ii) 
    229       ENDIF 
    230       IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1 
    231       ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj 
    232       ELSE                           ;   ij  = mj1(ij) 
    233       ENDIF 
    234       ! 
    235       IF( ii == jpi )   ii = ii-1       
    236       IF( ij == jpj )   ij = ij-1 
    237       ! 
    238       icb_utl_bilin_h = ( pfld(ii,ij  ) * (1.-zi) + pfld(ii+1,ij  ) * zi ) * (1.-zj)   & 
    239          &            + ( pfld(ii,ij+1) * (1.-zi) + pfld(ii+1,ij+1) * zi ) *     zj 
    240       ! 
    241    END FUNCTION icb_utl_bilin_h 
    242  
    243120 
    244121   REAL(wp) FUNCTION icb_utl_bilin( pfld, pi, pj, cd_type ) 
     
    270147            zj = pj - REAL(ij,wp) 
    271148         CASE ( 'U' ) 
    272             ii = MAX(1, INT( pi-0.5 )) 
    273             ij = MAX(1, INT( pj     ))    ! U-point 
    274             zi = pi - 0.5 - REAL(ii,wp) 
     149            ii = MAX(1, INT( pi-0.5_wp )) 
     150            ij = MAX(1, INT( pj        ))    ! U-point 
     151            zi = pi - 0.5_wp - REAL(ii,wp) 
    275152            zj = pj - REAL(ij,wp) 
    276153         CASE ( 'V' ) 
    277             ii = MAX(1, INT( pi     )) 
    278             ij = MAX(1, INT( pj-0.5 ))    ! V-point 
     154            ii = MAX(1, INT( pi        )) 
     155            ij = MAX(1, INT( pj-0.5_wp ))    ! V-point 
    279156            zi = pi - REAL(ii,wp) 
    280             zj = pj - 0.5 - REAL(ij,wp) 
     157            zj = pj - 0.5_wp - REAL(ij,wp) 
    281158         CASE ( 'F' ) 
    282             ii = MAX(1, INT( pi-0.5 )) 
    283             ij = MAX(1, INT( pj-0.5 ))    ! F-point 
    284             zi = pi - 0.5 - REAL(ii,wp) 
    285             zj = pj - 0.5 - REAL(ij,wp) 
     159            ii = MAX(1, INT( pi-0.5_wp )) 
     160            ij = MAX(1, INT( pj-0.5_wp ))    ! F-point 
     161            zi = pi - 0.5_wp - REAL(ii,wp) 
     162            zj = pj - 0.5_wp - REAL(ij,wp) 
    286163      END SELECT 
    287164      ! 
    288       ! find position in this processor. Prevent near edge problems (see #1389) 
    289       IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1 
    290       ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi 
    291       ELSE                           ;   ii = mi1(ii) 
    292       ENDIF 
    293       IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1 
    294       ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj 
    295       ELSE                           ;   ij  = mj1(ij) 
    296       ENDIF 
    297       ! 
    298       IF( ii == jpi )   ii = ii-1       
    299       IF( ij == jpj )   ij = ij-1 
     165      ! find position in this processor. 
     166      ii = mi1(ii) ; ij  = mj1(ij) 
    300167      ! 
    301168      icb_utl_bilin = ( pfld(ii,ij  ) * (1.-zi) + pfld(ii+1,ij  ) * zi ) * (1.-zj)   & 
     
    333200      zj = pj - REAL(ij,wp) 
    334201      ! 
    335       ! find position in this processor. Prevent near edge problems (see #1389) 
    336       IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1 
    337       ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi 
    338       ELSE                           ;   ii = mi1(ii) 
    339       ENDIF 
    340       IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1 
    341       ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj 
    342       ELSE                           ;   ij  = mj1(ij) 
    343       ENDIF 
    344       ! 
    345       IF( ii == jpi )   ii = ii-1       
    346       IF( ij == jpj )   ij = ij-1 
     202      ! find position in this processor. 
     203      ii = mi1(ii) ; ij  = mj1(ij) 
    347204      ! 
    348205      z4(1) = pfld(ii  ,ij  ) 
     
    392249      zi = pi - REAL(ii,wp)          !!gm use here mig, mjg arrays 
    393250      zj = pj - REAL(ij,wp) 
    394  
    395       ! find position in this processor. Prevent near edge problems (see #1389) 
    396       IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1 
    397       ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi 
    398       ELSE                           ;   ii = mi1(ii) 
    399       ENDIF 
    400       IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1 
    401       ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj 
    402       ELSE                           ;   ij  = mj1(ij) 
    403       ENDIF 
    404       ! 
    405       IF( ii == jpi )   ii = ii-1       
    406       IF( ij == jpj )   ij = ij-1 
     251      ! 
     252      ! find position in this processor. 
     253      ii = mi1(ii) ; ij  = mj1(ij) 
    407254      ! 
    408255      IF(    0.0_wp <= zi .AND. zi < 0.5_wp   ) THEN 
  • NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/LBC/lbclnk.F90

    r10425 r10679  
    4141   END INTERFACE 
    4242   ! 
    43    INTERFACE lbc_lnk_icb 
    44       MODULE PROCEDURE mpp_lnk_2d_icb 
    45    END INTERFACE 
    4643 
    4744   PUBLIC   lbc_lnk       ! ocean/ice lateral boundary conditions 
    4845   PUBLIC   lbc_lnk_multi ! modified ocean/ice lateral boundary conditions 
    4946   PUBLIC   lbc_bdy_lnk   ! ocean lateral BDY boundary conditions 
    50    PUBLIC   lbc_lnk_icb   ! iceberg lateral boundary conditions 
    5147 
    5248   !!---------------------------------------------------------------------- 
     
    9389   END INTERFACE 
    9490   ! 
    95    INTERFACE lbc_lnk_icb 
    96       MODULE PROCEDURE lbc_lnk_2d_icb 
    97    END INTERFACE 
    9891    
    9992   PUBLIC   lbc_lnk       ! ocean/ice  lateral boundary conditions 
    10093   PUBLIC   lbc_lnk_multi ! modified ocean/ice lateral boundary conditions 
    10194   PUBLIC   lbc_bdy_lnk   ! ocean lateral BDY boundary conditions 
    102    PUBLIC   lbc_lnk_icb   ! iceberg lateral boundary conditions 
    10395    
    10496   !!---------------------------------------------------------------------- 
     
    214206 
    215207 
    216 !!gm  This routine should be removed with an optional halos size added in argument of generic routines 
    217  
    218    SUBROUTINE lbc_lnk_2d_icb( cdname, pt2d, cd_type, psgn, ki, kj ) 
    219       !!---------------------------------------------------------------------- 
    220       CHARACTER(len=*)        , INTENT(in   ) ::   cdname      ! name of the calling subroutine 
    221       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt2d      ! 2D array on which the lbc is applied 
    222       CHARACTER(len=1)        , INTENT(in   ) ::   cd_type   ! nature of pt3d grid-points 
    223       REAL(wp)                , INTENT(in   ) ::   psgn      ! sign used across north fold  
    224       INTEGER                 , INTENT(in   ) ::   ki, kj    ! sizes of extra halo (not needed in non-mpp) 
    225       !!---------------------------------------------------------------------- 
    226       CALL lbc_lnk_2d( cdname, pt2d, cd_type, psgn ) 
    227    END SUBROUTINE lbc_lnk_2d_icb 
    228 !!gm end 
    229208 
    230209#endif 
  • NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/LBC/lib_mpp.F90

    r10538 r10679  
    4141   !!   mynode        : indentify the processor unit 
    4242   !!   mpp_lnk       : interface (defined in lbclnk) for message passing of 2d or 3d arrays (mpp_lnk_2d, mpp_lnk_3d) 
    43    !!   mpp_lnk_icb   : interface for message passing of 2d arrays with extra halo for icebergs (mpp_lnk_2d_icb) 
    4443   !!   mpprecv       : 
    4544   !!   mppsend       : 
     
    5453   !!   mppstop       : 
    5554   !!   mpp_ini_north : initialisation of north fold 
    56    !!   mpp_lbc_north_icb : alternative to mpp_nfd for extra outer halo with icebergs 
    5755   !!---------------------------------------------------------------------- 
    5856   USE dom_oce        ! ocean space and time domain 
     
    8078   PUBLIC   mynode, mppstop, mppsync, mpp_comm_free 
    8179   PUBLIC   mpp_ini_north 
    82    PUBLIC   mpp_lnk_2d_icb 
    83    PUBLIC   mpp_lbc_north_icb 
    8480   PUBLIC   mpp_min, mpp_max, mpp_sum, mpp_minloc, mpp_maxloc 
    8581   PUBLIC   mpp_delay_max, mpp_delay_sum, mpp_delay_rcv 
     
    11851181      ! 
    11861182   END SUBROUTINE DDPDD_MPI 
    1187  
    1188  
    1189    SUBROUTINE mpp_lbc_north_icb( pt2d, cd_type, psgn, kextj) 
    1190       !!--------------------------------------------------------------------- 
    1191       !!                   ***  routine mpp_lbc_north_icb  *** 
    1192       !! 
    1193       !! ** Purpose :   Ensure proper north fold horizontal bondary condition 
    1194       !!              in mpp configuration in case of jpn1 > 1 and for 2d 
    1195       !!              array with outer extra halo 
    1196       !! 
    1197       !! ** Method  :   North fold condition and mpp with more than one proc 
    1198       !!              in i-direction require a specific treatment. We gather 
    1199       !!              the 4+kextj northern lines of the global domain on 1 
    1200       !!              processor and apply lbc north-fold on this sub array. 
    1201       !!              Then we scatter the north fold array back to the processors. 
    1202       !!              This routine accounts for an extra halo with icebergs 
    1203       !!              and assumes ghost rows and columns have been suppressed. 
    1204       !! 
    1205       !!---------------------------------------------------------------------- 
    1206       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt2d     ! 2D array with extra halo 
    1207       CHARACTER(len=1)        , INTENT(in   ) ::   cd_type  ! nature of pt3d grid-points 
    1208       !                                                     !   = T ,  U , V , F or W -points 
    1209       REAL(wp)                , INTENT(in   ) ::   psgn     ! = -1. the sign change across the 
    1210       !!                                                    ! north fold, =  1. otherwise 
    1211       INTEGER                 , INTENT(in   ) ::   kextj    ! Extra halo width at north fold 
    1212       ! 
    1213       INTEGER ::   ji, jj, jr 
    1214       INTEGER ::   ierr, itaille, ildi, ilei, iilb 
    1215       INTEGER ::   ipj, ij, iproc 
    1216       ! 
    1217       REAL(wp), DIMENSION(:,:)  , ALLOCATABLE  ::  ztab_e, znorthloc_e 
    1218       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE  ::  znorthgloio_e 
    1219       !!---------------------------------------------------------------------- 
    1220       ! 
    1221       ipj=4 
    1222       ALLOCATE(        ztab_e(jpiglo, 1-kextj:ipj+kextj)       ,       & 
    1223      &            znorthloc_e(jpimax, 1-kextj:ipj+kextj)       ,       & 
    1224      &          znorthgloio_e(jpimax, 1-kextj:ipj+kextj,jpni)    ) 
    1225       ! 
    1226       ztab_e(:,:)      = 0._wp 
    1227       znorthloc_e(:,:) = 0._wp 
    1228       ! 
    1229       ij = 1 - kextj 
    1230       ! put the last ipj+2*kextj lines of pt2d into znorthloc_e  
    1231       DO jj = jpj - ipj + 1 - kextj , jpj + kextj 
    1232          znorthloc_e(1:jpi,ij)=pt2d(1:jpi,jj) 
    1233          ij = ij + 1 
    1234       END DO 
    1235       ! 
    1236       itaille = jpimax * ( ipj + 2*kextj ) 
    1237       ! 
    1238       IF( ln_timing ) CALL tic_tac(.TRUE.) 
    1239       CALL MPI_ALLGATHER( znorthloc_e(1,1-kextj)    , itaille, MPI_DOUBLE_PRECISION,    & 
    1240          &                znorthgloio_e(1,1-kextj,1), itaille, MPI_DOUBLE_PRECISION,    & 
    1241          &                ncomm_north, ierr ) 
    1242       ! 
    1243       IF( ln_timing ) CALL tic_tac(.FALSE.) 
    1244       ! 
    1245       DO jr = 1, ndim_rank_north            ! recover the global north array 
    1246          iproc = nrank_north(jr) + 1 
    1247          ildi = nldit (iproc) 
    1248          ilei = nleit (iproc) 
    1249          iilb = nimppt(iproc) 
    1250          DO jj = 1-kextj, ipj+kextj 
    1251             DO ji = ildi, ilei 
    1252                ztab_e(ji+iilb-1,jj) = znorthgloio_e(ji,jj,jr) 
    1253             END DO 
    1254          END DO 
    1255       END DO 
    1256  
    1257       ! 2. North-Fold boundary conditions 
    1258       ! ---------------------------------- 
    1259       CALL lbc_nfd( ztab_e(:,1-kextj:ipj+kextj), cd_type, psgn, kextj ) 
    1260  
    1261       ij = 1 - kextj 
    1262       !! Scatter back to pt2d 
    1263       DO jj = jpj - ipj + 1 - kextj , jpj + kextj 
    1264          DO ji= 1, jpi 
    1265             pt2d(ji,jj) = ztab_e(ji+nimpp-1,ij) 
    1266          END DO 
    1267          ij  = ij +1 
    1268       END DO 
    1269       ! 
    1270       DEALLOCATE( ztab_e, znorthloc_e, znorthgloio_e ) 
    1271       ! 
    1272    END SUBROUTINE mpp_lbc_north_icb 
    1273  
    1274  
    1275    SUBROUTINE mpp_lnk_2d_icb( cdname, pt2d, cd_type, psgn, kexti, kextj ) 
    1276       !!---------------------------------------------------------------------- 
    1277       !!                  ***  routine mpp_lnk_2d_icb  *** 
    1278       !! 
    1279       !! ** Purpose :   Message passing management for 2d array (with extra halo for icebergs) 
    1280       !!                This routine receives a (1-kexti:jpi+kexti,1-kexti:jpj+kextj) 
    1281       !!                array (usually (0:jpi+1, 0:jpj+1)) from lbc_lnk_icb calls. 
    1282       !! 
    1283       !! ** Method  :   Use mppsend and mpprecv function for passing mask 
    1284       !!      between processors following neighboring subdomains. 
    1285       !!            domain parameters 
    1286       !!                    jpi    : first dimension of the local subdomain 
    1287       !!                    jpj    : second dimension of the local subdomain 
    1288       !!                    kexti  : number of columns for extra outer halo 
    1289       !!                    kextj  : number of rows for extra outer halo 
    1290       !!                    nbondi : mark for "east-west local boundary" 
    1291       !!                    nbondj : mark for "north-south local boundary" 
    1292       !!                    noea   : number for local neighboring processors 
    1293       !!                    nowe   : number for local neighboring processors 
    1294       !!                    noso   : number for local neighboring processors 
    1295       !!                    nono   : number for local neighboring processors 
    1296       !!---------------------------------------------------------------------- 
    1297       CHARACTER(len=*)                                        , INTENT(in   ) ::   cdname      ! name of the calling subroutine 
    1298       REAL(wp), DIMENSION(1-kexti:jpi+kexti,1-kextj:jpj+kextj), INTENT(inout) ::   pt2d     ! 2D array with extra halo 
    1299       CHARACTER(len=1)                                        , INTENT(in   ) ::   cd_type  ! nature of ptab array grid-points 
    1300       REAL(wp)                                                , INTENT(in   ) ::   psgn     ! sign used across the north fold 
    1301       INTEGER                                                 , INTENT(in   ) ::   kexti    ! extra i-halo width 
    1302       INTEGER                                                 , INTENT(in   ) ::   kextj    ! extra j-halo width 
    1303       ! 
    1304       INTEGER  ::   jl   ! dummy loop indices 
    1305       INTEGER  ::   imigr, iihom, ijhom        ! local integers 
    1306       INTEGER  ::   ipreci, iprecj             !   -       - 
    1307       INTEGER  ::   ml_req1, ml_req2, ml_err   ! for key_mpi_isend 
    1308       INTEGER, DIMENSION(MPI_STATUS_SIZE) ::   ml_stat   ! for key_mpi_isend 
    1309       !! 
    1310       REAL(wp), DIMENSION(1-kexti:jpi+kexti,nn_hls+kextj,2) ::   r2dns, r2dsn 
    1311       REAL(wp), DIMENSION(1-kextj:jpj+kextj,nn_hls+kexti,2) ::   r2dwe, r2dew 
    1312       !!---------------------------------------------------------------------- 
    1313  
    1314       ipreci = nn_hls + kexti      ! take into account outer extra 2D overlap area 
    1315       iprecj = nn_hls + kextj 
    1316  
    1317       IF( narea == 1 .AND. numcom == -1 ) CALL mpp_report( cdname, 1, 1, 1, ld_lbc = .TRUE. ) 
    1318  
    1319       ! 1. standard boundary treatment 
    1320       ! ------------------------------ 
    1321       ! Order matters Here !!!! 
    1322       ! 
    1323       !                                      ! East-West boundaries 
    1324       !                                           !* Cyclic east-west 
    1325       IF( l_Iperio ) THEN 
    1326          pt2d(1-kexti:     1   ,:) = pt2d(jpim1-kexti: jpim1 ,:)       ! east 
    1327          pt2d(  jpi  :jpi+kexti,:) = pt2d(     2     :2+kexti,:)       ! west 
    1328          ! 
    1329       ELSE                                        !* closed 
    1330          IF( .NOT. cd_type == 'F' )   pt2d(  1-kexti   :nn_hls   ,:) = 0._wp    ! east except at F-point 
    1331                                       pt2d(jpi-nn_hls+1:jpi+kexti,:) = 0._wp    ! west 
    1332       ENDIF 
    1333       !                                      ! North-South boundaries 
    1334       IF( l_Jperio ) THEN                         !* cyclic (only with no mpp j-split) 
    1335          pt2d(:,1-kextj:     1   ) = pt2d(:,jpjm1-kextj:  jpjm1)       ! north 
    1336          pt2d(:,  jpj  :jpj+kextj) = pt2d(:,     2     :2+kextj)       ! south 
    1337       ELSE                                        !* closed 
    1338          IF( .NOT. cd_type == 'F' )   pt2d(:,  1-kextj   :nn_hls   ) = 0._wp    ! north except at F-point 
    1339                                       pt2d(:,jpj-nn_hls+1:jpj+kextj) = 0._wp    ! south 
    1340       ENDIF 
    1341       ! 
    1342  
    1343       ! north fold treatment 
    1344       ! ----------------------- 
    1345       IF( npolj /= 0 ) THEN 
    1346          ! 
    1347          SELECT CASE ( jpni ) 
    1348                    CASE ( 1 )     ;   CALL lbc_nfd          ( pt2d(1:jpi,1:jpj+kextj), cd_type, psgn, kextj ) 
    1349                    CASE DEFAULT   ;   CALL mpp_lbc_north_icb( pt2d(1:jpi,1:jpj+kextj), cd_type, psgn, kextj ) 
    1350          END SELECT 
    1351          ! 
    1352       ENDIF 
    1353  
    1354       ! 2. East and west directions exchange 
    1355       ! ------------------------------------ 
    1356       ! we play with the neigbours AND the row number because of the periodicity 
    1357       ! 
    1358       SELECT CASE ( nbondi )      ! Read Dirichlet lateral conditions 
    1359       CASE ( -1, 0, 1 )                ! all exept 2 (i.e. close case) 
    1360          iihom = jpi-nreci-kexti 
    1361          DO jl = 1, ipreci 
    1362             r2dew(:,jl,1) = pt2d(nn_hls+jl,:) 
    1363             r2dwe(:,jl,1) = pt2d(iihom +jl,:) 
    1364          END DO 
    1365       END SELECT 
    1366       ! 
    1367       !                           ! Migrations 
    1368       imigr = ipreci * ( jpj + 2*kextj ) 
    1369       ! 
    1370       IF( ln_timing ) CALL tic_tac(.TRUE.) 
    1371       ! 
    1372       SELECT CASE ( nbondi ) 
    1373       CASE ( -1 ) 
    1374          CALL mppsend( 2, r2dwe(1-kextj,1,1), imigr, noea, ml_req1 ) 
    1375          CALL mpprecv( 1, r2dew(1-kextj,1,2), imigr, noea ) 
    1376          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 
    1377       CASE ( 0 ) 
    1378          CALL mppsend( 1, r2dew(1-kextj,1,1), imigr, nowe, ml_req1 ) 
    1379          CALL mppsend( 2, r2dwe(1-kextj,1,1), imigr, noea, ml_req2 ) 
    1380          CALL mpprecv( 1, r2dew(1-kextj,1,2), imigr, noea ) 
    1381          CALL mpprecv( 2, r2dwe(1-kextj,1,2), imigr, nowe ) 
    1382          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 
    1383          IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err) 
    1384       CASE ( 1 ) 
    1385          CALL mppsend( 1, r2dew(1-kextj,1,1), imigr, nowe, ml_req1 ) 
    1386          CALL mpprecv( 2, r2dwe(1-kextj,1,2), imigr, nowe ) 
    1387          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 
    1388       END SELECT 
    1389       ! 
    1390       IF( ln_timing ) CALL tic_tac(.FALSE.) 
    1391       ! 
    1392       !                           ! Write Dirichlet lateral conditions 
    1393       iihom = jpi - nn_hls 
    1394       ! 
    1395       SELECT CASE ( nbondi ) 
    1396       CASE ( -1 ) 
    1397          DO jl = 1, ipreci 
    1398             pt2d(iihom+jl,:) = r2dew(:,jl,2) 
    1399          END DO 
    1400       CASE ( 0 ) 
    1401          DO jl = 1, ipreci 
    1402             pt2d(jl-kexti,:) = r2dwe(:,jl,2) 
    1403             pt2d(iihom+jl,:) = r2dew(:,jl,2) 
    1404          END DO 
    1405       CASE ( 1 ) 
    1406          DO jl = 1, ipreci 
    1407             pt2d(jl-kexti,:) = r2dwe(:,jl,2) 
    1408          END DO 
    1409       END SELECT 
    1410  
    1411  
    1412       ! 3. North and south directions 
    1413       ! ----------------------------- 
    1414       ! always closed : we play only with the neigbours 
    1415       ! 
    1416       IF( nbondj /= 2 ) THEN      ! Read Dirichlet lateral conditions 
    1417          ijhom = jpj-nrecj-kextj 
    1418          DO jl = 1, iprecj 
    1419             r2dsn(:,jl,1) = pt2d(:,ijhom +jl) 
    1420             r2dns(:,jl,1) = pt2d(:,nn_hls+jl) 
    1421          END DO 
    1422       ENDIF 
    1423       ! 
    1424       !                           ! Migrations 
    1425       imigr = iprecj * ( jpi + 2*kexti ) 
    1426       ! 
    1427       IF( ln_timing ) CALL tic_tac(.TRUE.) 
    1428       ! 
    1429       SELECT CASE ( nbondj ) 
    1430       CASE ( -1 ) 
    1431          CALL mppsend( 4, r2dsn(1-kexti,1,1), imigr, nono, ml_req1 ) 
    1432          CALL mpprecv( 3, r2dns(1-kexti,1,2), imigr, nono ) 
    1433          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 
    1434       CASE ( 0 ) 
    1435          CALL mppsend( 3, r2dns(1-kexti,1,1), imigr, noso, ml_req1 ) 
    1436          CALL mppsend( 4, r2dsn(1-kexti,1,1), imigr, nono, ml_req2 ) 
    1437          CALL mpprecv( 3, r2dns(1-kexti,1,2), imigr, nono ) 
    1438          CALL mpprecv( 4, r2dsn(1-kexti,1,2), imigr, noso ) 
    1439          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 
    1440          IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err) 
    1441       CASE ( 1 ) 
    1442          CALL mppsend( 3, r2dns(1-kexti,1,1), imigr, noso, ml_req1 ) 
    1443          CALL mpprecv( 4, r2dsn(1-kexti,1,2), imigr, noso ) 
    1444          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 
    1445       END SELECT 
    1446       ! 
    1447       IF( ln_timing ) CALL tic_tac(.FALSE.) 
    1448       ! 
    1449       !                           ! Write Dirichlet lateral conditions 
    1450       ijhom = jpj - nn_hls 
    1451       ! 
    1452       SELECT CASE ( nbondj ) 
    1453       CASE ( -1 ) 
    1454          DO jl = 1, iprecj 
    1455             pt2d(:,ijhom+jl) = r2dns(:,jl,2) 
    1456          END DO 
    1457       CASE ( 0 ) 
    1458          DO jl = 1, iprecj 
    1459             pt2d(:,jl-kextj) = r2dsn(:,jl,2) 
    1460             pt2d(:,ijhom+jl) = r2dns(:,jl,2) 
    1461          END DO 
    1462       CASE ( 1 ) 
    1463          DO jl = 1, iprecj 
    1464             pt2d(:,jl-kextj) = r2dsn(:,jl,2) 
    1465          END DO 
    1466       END SELECT 
    1467       ! 
    1468    END SUBROUTINE mpp_lnk_2d_icb 
    1469  
    14701183 
    14711184   SUBROUTINE mpp_report( cdname, kpk, kpl, kpf, ld_lbc, ld_glb, ld_dlg ) 
Note: See TracChangeset for help on using the changeset viewer.