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 10679 for NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB – NEMO

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

branch for solution 2 of ticket #2238

Location:
NEMO/branches/2019/fix_ticket2238_solution2
Files:
5 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 
Note: See TracChangeset for help on using the changeset viewer.