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 14573 – NEMO

Changeset 14573


Ignore:
Timestamp:
2021-03-03T14:39:55+01:00 (3 years ago)
Author:
antsia
Message:

Merge with dev_isf_remapping_UKESM_GO6package_r9314@11248

Location:
branches/NERC/dev_r5518_GO6_r9321_isf_merge_UKESM_AIS/NEMOGCM
Files:
5 edited
1 copied

Legend:

Unmodified
Added
Removed

  • branches/NERC/dev_r5518_GO6_r9321_isf_merge_UKESM_AIS/NEMOGCM/CONFIG/SHARED/namelist_ref

    r14572 r14573  
    116116                        !!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above] 
    117117   rn_thetb    =    1.0    !  bottom control parameter  (0<=thetb<= 1) 
     118/ 
     119!----------------------------------------------------------------------- 
     120&namzgr_isf    !   isf cavity geometry definition 
     121!----------------------------------------------------------------------- 
     122   rn_isfdep_min    = 10.         ! minimum isf draft tickness (if lower, isf draft set to this value) 
     123   rn_glhw_min      = 1.e-3       ! minimum water column thickness to define the grounding line 
     124   rn_isfhw_min     = 10          ! minimum water column thickness in the cavity once the grounding line defined. 
     125   ln_isfchannel    = .false.     ! remove channel (based on 2d mask build from isfdraft-bathy) 
     126   ln_isfconnect    = .false.     ! force connection under the ice shelf (based on 2d mask build from isfdraft-bathy) 
     127   nn_kisfmax       = 999         ! limiter in level on the previous condition. (if change larger than this number, get back to value before we enforce the connection) 
     128   rn_zisfmax       = 7000.       ! limiter in m     on the previous condition. (if change larger than this number, get back to value before we enforce the connection) 
     129   ln_isfcheminey   = .false.     ! close cheminey 
     130   ln_isfsubgl      = .true.      ! remove subglacial lake created by the remapping process 
     131   rn_isfsubgllon   =    0.0      !  longitude of the seed to determine the open ocean 
     132   rn_isfsubgllat   =    0.0      !  latitude  of the seed to determine the open ocean 
    118133/ 
    119134!----------------------------------------------------------------------- 

  • branches/NERC/dev_r5518_GO6_r9321_isf_merge_UKESM_AIS/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

    r14572 r14573  
    2525   USE oce             ! ocean dynamics and tracers 
    2626   USE dom_oce         ! ocean space and time domain 
     27   USE domngb          ! find nearest wet point 
     28   USE domutl          ! fill closed 3d pool below isf 
     29   USE domzgr, ONLY : ln_isfsubgl, rn_isfsubgllon, rn_isfsubgllat ! import ln_isfsubgl to mask close sea below isf 
     30   ! 
    2731   USE in_out_manager  ! I/O manager 
    2832   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    … …  
    133137      INTEGER  ::   iif, iil, ii0, ii1, ii   ! local integers 
    134138      INTEGER  ::   ijf, ijl, ij0, ij1       !   -       - 
     139      INTEGER  ::   jiseed, jjseed           !   -       - 
    135140      INTEGER  ::   ios 
    136141      INTEGER  ::   isrow                    ! index for ORCA1 starting row 
    … …  
    187192      END DO   
    188193       
    189       ! (ISF) define barotropic mask and mask the ice shelf point 
    190       ssmask(:,:)=tmask(:,:,1) ! at this stage ice shelf is not masked 
    191        
    192194      DO jk = 1, jpk 
    193195         DO jj = 1, jpj 
    … …  
    199201         END DO   
    200202      END DO   
     203      ! 
     204      ! (ISF) define barotropic mask and mask the ice shelf point 
     205      DO jj = 1, jpj 
     206         DO ji = 1, jpi   ! vector loop 
     207            ssmask(ji,jj)  = MIN(1._wp,SUM(tmask(ji,jj,:))) 
     208         END DO 
     209      END DO 
     210      ! 
     211      IF ( ln_isfsubgl ) THEN 
     212         ! check closed wet pool 
     213         CALL dom_ngb(rn_isfsubgllon, rn_isfsubgllat, jiseed, jjseed, 'T', lwet=.TRUE.) 
     214         CALL fill_pool( jiseed, jjseed, 1, tmask, -1._wp ) 
     215         ! at this point itab3d (:,1:ijmax,:) can have 3 different values : 
     216         !              0 where there where already 0 
     217         !              -1 where the ocean points are connected 
     218         !              1 where ocean points in tmask are not connected 
     219         IF (lwp) THEN 
     220            WRITE(numout,*) 
     221            WRITE(numout,*)'dommsk : removal of subglacial lakes ' 
     222            WRITE(numout,*)'~~~~~~~' 
     223            WRITE(numout,*)'Number of disconected points : ', COUNT(  (tmask(:,:,:) == 1) ) 
     224            WRITE(numout,*)'lon/lat seed to detect main ocean is: ', rn_isfsubgllon, rn_isfsubgllat 
     225            WRITE(numout,*)'i/j     seed to detect main ocean is: ', jiseed, jjseed 
     226         END IF 
     227         DO jk = 1, jpk 
     228            WHERE (tmask(:,:,jk) > 0 .AND. misfdep(:,:) > 1) tmask(:,:,jk) = 0 ! remove only subglacial lake (ie similar to close sea only below an ice shelf  
     229            WHERE (tmask(:,:,jk) < 0) tmask(:,:,jk) = 1 ! restore mask value 
     230         END DO 
     231         ! update ssmask 
     232         DO jj = 1, jpj 
     233            DO ji = 1, jpi   ! vector loop 
     234               ssmask(ji,jj)  = MIN(1._wp,SUM(tmask(ji,jj,:))) 
     235            END DO 
     236         END DO 
     237         ! update mbathy, misfdep, bathy, risfdep 
     238         bathy(:,:)   = bathy(:,:)   * ssmask(:,:) 
     239         risfdep(:,:) = risfdep(:,:) * ssmask(:,:) 
     240         WHERE ( ssmask(:,:) == 0._wp ) 
     241            misfdep(:,:) = 1 
     242            mbathy(:,:)  = 0 
     243         END WHERE 
     244      END IF 
    201245 
    202246!!gm  ???? 

  • branches/NERC/dev_r5518_GO6_r9321_isf_merge_UKESM_AIS/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90

    r6486 r14573  
    2828CONTAINS 
    2929 
    30    SUBROUTINE dom_ngb( plon, plat, kii, kjj, cdgrid ) 
     30   SUBROUTINE dom_ngb( plon, plat, kii, kjj, cdgrid, lwet, kkk ) 
    3131      !!---------------------------------------------------------------------- 
    3232      !!                    ***  ROUTINE dom_ngb  *** 
    3333      !! 
    34       !! ** Purpose :   find the closest grid point from a given lon/lat position 
     34      !! ** Purpose :   find the closest wet grid point from a given lon/lat position 
    3535      !! 
    3636      !! ** Method  :   look for minimum distance in cylindrical projection  
    … …  
    4040      REAL(wp)        , INTENT(in   ) ::   plon, plat   ! longitude,latitude of the point 
    4141      INTEGER         , INTENT(  out) ::   kii, kjj     ! i-,j-index of the closes grid point 
     42      INTEGER         , INTENT(in   ), OPTIONAL :: kkk  ! k-index of the mask level used 
     43      LOGICAL         , INTENT(in   ), OPTIONAL :: lwet ! logical to decide if we look for every where (false)  
     44                                                                !    or only over the wet point (true) 
    4245      CHARACTER(len=1), INTENT(in   ) ::   cdgrid       ! grid name 'T', 'U', 'V', 'W' 
    4346      ! 
     47      INTEGER :: ik         ! working level 
    4448      INTEGER , DIMENSION(2) ::   iloc 
    4549      REAL(wp)               ::   zlon, zmini 
    4650      REAL(wp), POINTER, DIMENSION(:,:) ::  zglam, zgphi, zmask, zdist 
     51      LOGICAL :: llwet      ! working logical 
    4752      !!-------------------------------------------------------------------- 
    4853      ! 
    … …  
    5156      CALL wrk_alloc( jpi, jpj, zglam, zgphi, zmask, zdist ) 
    5257      ! 
     58      ! select lat/lon variable 
     59      SELECT CASE( cdgrid ) 
     60         CASE( 'U' )  ; zglam(:,:) = glamu(:,:) ; zgphi(:,:) = gphiu(:,:) 
     61         CASE( 'V' )  ; zglam(:,:) = glamv(:,:) ; zgphi(:,:) = gphiv(:,:) 
     62         CASE( 'F' )  ; zglam(:,:) = glamf(:,:) ; zgphi(:,:) = gphif(:,:) 
     63         CASE DEFAULT ; zglam(:,:) = glamt(:,:) ; zgphi(:,:) = gphit(:,:) 
     64      END SELECT 
     65      ! 
     66      ! select vertical level 
     67      ik   = 1 
     68      IF ( PRESENT(kkk) ) ik=kkk 
     69      ! 
     70      ! select mask variable 
    5371      zmask(:,:) = 0._wp 
    54       SELECT CASE( cdgrid ) 
    55       CASE( 'U' )  ; zglam(:,:) = glamu(:,:) ; zgphi(:,:) = gphiu(:,:) ; zmask(nldi:nlei,nldj:nlej) = umask(nldi:nlei,nldj:nlej,1) 
    56       CASE( 'V' )  ; zglam(:,:) = glamv(:,:) ; zgphi(:,:) = gphiv(:,:) ; zmask(nldi:nlei,nldj:nlej) = vmask(nldi:nlei,nldj:nlej,1) 
    57       CASE( 'F' )  ; zglam(:,:) = glamf(:,:) ; zgphi(:,:) = gphif(:,:) ; zmask(nldi:nlei,nldj:nlej) = fmask(nldi:nlei,nldj:nlej,1) 
    58       CASE DEFAULT ; zglam(:,:) = glamt(:,:) ; zgphi(:,:) = gphit(:,:) ; zmask(nldi:nlei,nldj:nlej) = tmask(nldi:nlei,nldj:nlej,1) 
    59       END SELECT 
    60  
    61       zlon       = MOD( plon       + 720., 360. )                                     ! plon between    0 and 360 
    62       zglam(:,:) = MOD( zglam(:,:) + 720., 360. )                                     ! glam between    0 and 360 
    63       IF( zlon > 270. )   zlon = zlon - 360.                                          ! zlon between  -90 and 270 
    64       IF( zlon <  90. )   WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360.   ! glam between -180 and 180 
    65  
    66       zglam(:,:) = zglam(:,:) - zlon 
     72      llwet = .TRUE. 
     73      IF ( PRESENT(lwet)) llwet=lwet  
     74      IF (llwet) THEN 
     75         SELECT CASE( cdgrid ) 
     76            CASE( 'U' )  ; zmask(nldi:nlei,nldj:nlej) = umask(nldi:nlei,nldj:nlej,ik) 
     77            CASE( 'V' )  ; zmask(nldi:nlei,nldj:nlej) = vmask(nldi:nlei,nldj:nlej,ik) 
     78            CASE( 'F' )  ; zmask(nldi:nlei,nldj:nlej) = fmask(nldi:nlei,nldj:nlej,ik) 
     79            CASE DEFAULT ; zmask(nldi:nlei,nldj:nlej) = tmask(nldi:nlei,nldj:nlej,ik) 
     80         END SELECT 
     81      ELSE 
     82         zmask(nldi:nlei,nldj:nlej)=1.e0 
     83      END IF 
     84      ! 
     85      ! compute distance 
     86      IF (jphgr_msh /= 2 .AND. jphgr_msh /= 3) THEN 
     87         zlon       = MOD( plon       + 720., 360. )                                     ! plon between    0 and 360 
     88         zglam(:,:) = MOD( zglam(:,:) + 720., 360. )                                     ! glam between    0 and 360 
     89         IF( zlon > 270. )   zlon = zlon - 360.                                          ! zlon between  -90 and 270 
     90         IF( zlon <  90. )   WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360.   ! glam between -180 and 180 
     91         zglam(:,:) = zglam(:,:) - zlon 
     92      ELSE 
     93         zglam(:,:) = zglam(:,:) - plon 
     94      END IF 
     95      ! 
    6796      zgphi(:,:) = zgphi(:,:) - plat 
    6897      zdist(:,:) = zglam(:,:) * zglam(:,:) + zgphi(:,:) * zgphi(:,:) 
    69        
     98      ! 
     99      ! find min 
    70100      IF( lk_mpp ) THEN   
    71101         CALL mpp_minloc( zdist(:,:), zmask, zmini, kii, kjj) 

  • branches/NERC/dev_r5518_GO6_r9321_isf_merge_UKESM_AIS/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90

    r14572 r14573  
    205205      zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp ) 
    206206      CALL iom_rstput( 0, 0, inum4, 'isfdraft', zprt, ktype = jp_r4 )       !    ! nb of ocean T-points 
     207      zprt(:,:) = ssmask(:,:) * REAL( bathy(:,:) , wp ) 
     208      CALL iom_rstput( 0, 0, inum4, 'bathy', zprt, ktype = jp_r4 )       !    ! nb of ocean T-points 
     209             
    207210             
    208211      IF( ln_sco ) THEN                                         ! s-coordinate 

  • branches/NERC/dev_r5518_GO6_r9321_isf_merge_UKESM_AIS/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r6487 r14573  
    5050   PUBLIC   dom_zgr        ! called by dom_init.F90 
    5151 
     52   ! remove subglacial lake under the ice sheet. 
     53   LOGICAL , PUBLIC  :: ln_isfsubgl 
     54   REAL(wp), PUBLIC  :: rn_isfsubgllon, rn_isfsubgllat 
    5255   !                              !!* Namelist namzgr_sco * 
    5356   LOGICAL  ::   ln_s_sh94         ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T) 
    … …  
    395398      IF(lwp) WRITE(numout,*) '    zgr_bat : defines level and meter bathymetry' 
    396399      IF(lwp) WRITE(numout,*) '    ~~~~~~~' 
     400      ! 
     401      ! (ISF) initialisation ice shelf draft and top level 
     402      risfdep(:,:)=0._wp 
     403      misfdep(:,:)=1 
    397404      !                                               ! ================== !  
    398405      IF( ntopo == 0 .OR. ntopo == -1 ) THEN          !   defined by hand  ! 
    … …  
    484491            END DO 
    485492         END DO 
    486          risfdep(:,:)=0.e0 
    487          misfdep(:,:)=1 
    488493         ! 
    489494         DEALLOCATE( idta, zdta ) 
    … …  
    535540            CALL iom_close( inum ) 
    536541            !                                                 
    537             risfdep(:,:)=0._wp          
    538             misfdep(:,:)=1              
    539542            IF ( ln_isfcav ) THEN 
    540543               CALL iom_open ( 'isf_draft_meter.nc', inum )  
    … …  
    12601263      !! ** Purpose :   check the bathymetry in levels 
    12611264      !!    
    1262       !! ** Method  :   THe water column have to contained at least 2 cells 
     1265      !! ** Method  :   The water column have to contained at least 2 cells 
    12631266      !!                Bathymetry and isfdraft are modified (dig/close) to respect 
    12641267      !!                this criterion. 
    1265       !!                  
     1268      !!                          0.0 = definition of minal ice shelf draft 
     1269      !!                digging and filling base on depth criterion only 
     1270      !!                          1.0 = set iceshelf to the minimum depth allowed 
     1271      !!                          1.1 = ground ice shelf if water column less than X m 
     1272      !!                          1.2 = ensure a minimum thickness for iceshelf cavity 
     1273      !!                          1.3 = remove channels and single point 'bay' 
     1274      !!                          1.4 = close single isolated point 
     1275      !!                compute level 
     1276      !!                          2.0 = compute level 
     1277      !!                          2.1 = ensure misfdep is not below bathymetry after step 2.0 
     1278      !!                digging and filling base on level criterion only                 
     1279      !!                          3.0 = dig to fit the 2 water level criterion (closed pool possible after this step) 
     1280      !!                          3.1 = dig enough to ensure connectivity of all the cell beneath an ice shelf  
     1281      !!                                (most of the close pool should be remove after this step) 
     1282      !!                          3.2 = fill chimney 
    12661283      !!    
    12671284      !! ** Action  : - test compatibility between isfdraft and bathy  
    12681285      !!              - bathy and isfdraft are modified 
    12691286      !!---------------------------------------------------------------------- 
    1270       !!    
    1271       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    1272       INTEGER  ::   ik, it           ! temporary integers 
    1273       INTEGER  ::   id, jd, nprocd 
    1274       INTEGER  ::   icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1   ! (ISF) 
    1275       LOGICAL  ::   ll_print         ! Allow  control print for debugging 
    1276       REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
    1277       REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t 
    1278       REAL(wp) ::   zmax, zmin       ! Maximum and minimum depth 
    1279       REAL(wp) ::   zdiff            ! temporary scalar 
    1280       REAL(wp) ::   zrefdep          ! temporary scalar 
    1281       REAL(wp) ::   zbathydiff, zrisfdepdiff  ! isf temporary scalar 
    1282       REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 2D workspace (ISH) 
    1283       INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep         ! 2D workspace (ISH) 
     1287      !!   
     1288      INTEGER  ::   ji, jj, jk, jn 
     1289      INTEGER  ::   ios, icompt 
     1290      INTEGER  ::   ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1   ! (ISF) 
     1291      INTEGER  ::   zmbathyij, zmbathyip1, zmbathyim1, zmbathyjp1, zmbathyjm1 
     1292      INTEGER , POINTER, DIMENSION(:,:)   ::   zmisfdep, zmbathy         ! 2D workspace (ISH) 
     1293      ! 
     1294      REAL(wp) ::   zdepth, vmskjp1, vmskjm1, umskip1, umskim1, umskip1_r, umskim1_r, vmskjp1_r, vmskjm1_r 
     1295      REAL(wp) :: zisfdep_min 
     1296      REAL(wp), POINTER, DIMENSION(:,:)   ::   zdummy, zmask, zrisfdep   ! 2D workspace (ISH) 
     1297      ! 
     1298      !!----------------- 
     1299      ! NAMELIST 
     1300      INTEGER  :: nn_kisfmax    = 999.                              !   
     1301      REAL(wp) :: rn_isfdep_min = 10.0_wp                           ! ice shelf minimal thickness  
     1302      REAL(wp) :: rn_glhw_min   = 1.0e-3                             ! threshold on hw to define grounding line water 
     1303      REAL(wp) :: rn_isfhw_min  = 1.0e-3                            ! threshold on hw to define isf draft into the cavity 
     1304      REAL(wp) :: rn_isfshallow = 0.0_wp                            ! threshold to define shallow ice shelf cavity 
     1305      REAL(wp) :: rn_zisfmax    = 6000.0_wp                         ! maximun meter of ice we are allowed to dig to assure connectivity 
     1306      LOGICAL  :: ln_isfcheminey= .FALSE. , ln_isfconnect= .FALSE. , ln_isfchannel= .FALSE. !  remove cheminey, assure connectivity, remove channel in water column (on z space not level space) 
    12841307      !!--------------------------------------------------------------------- 
     1308      NAMELIST/namzgr_isf/nn_kisfmax, rn_zisfmax, & 
     1309              &           rn_isfdep_min, rn_glhw_min, rn_isfhw_min, & 
     1310              &           ln_isfcheminey, ln_isfconnect, ln_isfchannel, ln_isfsubgl, rn_isfsubgllon, rn_isfsubgllat 
    12851311      ! 
    12861312      IF( nn_timing == 1 )  CALL timing_start('zgr_isf') 
    12871313      ! 
    1288       CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep) 
     1314      CALL wrk_alloc( jpi, jpj, zdummy, zmask, zrisfdep) 
    12891315      CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy ) 
    1290  
    1291  
    1292       ! (ISF) compute misfdep 
    1293       WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1   
    1294       ELSEWHERE                      ;                          misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level  
    1295       END WHERE   
    1296  
    1297       ! Compute misfdep for ocean points (i.e. first wet level)  
     1316      ! 
     1317      ! 
     1318      REWIND( numnam_ref )              ! Namelist namzgr_isf in reference namelist : ice shelf geometry definition 
     1319      READ  ( numnam_ref, namzgr_isf, IOSTAT = ios, ERR = 901) 
     1320901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_isf in reference namelist', lwp ) 
     1321 
     1322      REWIND( numnam_cfg )              ! Namelist namzgr_sco in configuration namelist : ice shelf geometry definition 
     1323      READ  ( numnam_cfg, namzgr_isf, IOSTAT = ios, ERR = 902 ) 
     1324902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_isf in configuration namelist', lwp ) 
     1325      IF(lwm) WRITE ( numond, namzgr_isf ) 
     1326      ! 
     1327      IF (lwp) THEN 
     1328         WRITE(numout,*) ' nn_kisfmax       = ',nn_kisfmax           ! 
     1329         WRITE(numout,*) ' rn_zisfmax       = ',rn_zisfmax           ! 
     1330         WRITE(numout,*) ' rn_isfdep_min    = ',rn_isfdep_min        ! 
     1331         WRITE(numout,*) ' rn_isfhw_min     = ',rn_isfhw_min         ! 
     1332         WRITE(numout,*) ' rn_glhw_min      = ',rn_glhw_min          ! 
     1333         WRITE(numout,*) ' ln_isfcheminey   = ',ln_isfcheminey       ! 
     1334         WRITE(numout,*) ' ln_isfconnect    = ',ln_isfconnect        ! 
     1335         WRITE(numout,*) ' ln_isfchannel    = ',ln_isfchannel        ! 
     1336      END IF 
     1337      ! 
     1338      ! 0.0 variable definition of minimal ice shelf draft allowed 
     1339      ! an ice shelf draft have to be larger than the first level thickness 
     1340      zisfdep_min = MAX(rn_isfdep_min,e3t_1d(1)) 
     1341 
     1342      ! 1.0 set iceshelf to the minimum depth allowed 
     1343      WHERE(risfdep(:,:) > 0.0_wp .AND. risfdep(:,:) < zisfdep_min) 
     1344         risfdep(:,:)=zisfdep_min 
     1345      END WHERE 
     1346      !   
     1347      ! 1.1 ground ice shelf if water column less than X m => set the grounding 
     1348      ! line position 
     1349      WHERE( bathy(:,:) - risfdep(:,:) < rn_glhw_min )  
     1350         risfdep(:,:)=0.0 ; misfdep(:,:) = 1 
     1351         bathy  (:,:)=0.0 ; mbathy (:,:) = 0 
     1352      END WHERE 
     1353      ! 
     1354      ! 1.2 ensure a minimum thickness for iceshelf cavity => avoid to negative 
     1355      ! e3t if ssh + sum(e3t*tmask) < 0 
     1356      DO jj = 1, jpj 
     1357         DO ji = 1, jpi 
     1358            IF (bathy(ji,jj)-risfdep(ji,jj) < rn_isfhw_min) THEN 
     1359               risfdep(ji,jj) = bathy(ji,jj) - rn_isfhw_min 
     1360               ! sanity check on risfdep (if < zisfdep_min) => we ground it  
     1361               IF (risfdep(ji,jj) < zisfdep_min) THEN  
     1362                  risfdep(ji,jj)=0.0 ; misfdep(ji,jj) = 1 
     1363                  bathy  (ji,jj)=0.0 ; mbathy (ji,jj) = 0 
     1364               END IF 
     1365            END IF 
     1366         END DO 
     1367      END DO 
     1368      ! 
     1369      ! 1.3 Remove channels and single point 'bay'. 
     1370      ! channel could be created if connectivity is enforced later. 
     1371      IF (ln_isfchannel) THEN 
     1372         zmask(:,:)=1._wp 
     1373         WHERE (mbathy(:,:)==0.0) 
     1374            zmask(:,:)=0._wp 
     1375         END WHERE 
     1376         DO jj = 2, jpjm1 
     1377            DO ji = 2, jpim1 
     1378               IF(  (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN 
     1379                  umskip1=zmask(ji,jj)*zmask(ji+1,jj  )  ! 1 = connexion 
     1380                  vmskjp1=zmask(ji,jj)*zmask(ji  ,jj+1)  ! 1 = connexion 
     1381                  umskim1=zmask(ji,jj)*zmask(ji-1,jj  )  ! 1 = connexion 
     1382                  vmskjm1=zmask(ji,jj)*zmask(ji  ,jj-1)  ! 1 = connexion 
     1383               ! zonal channel and single bay 
     1384                  IF ((umskip1+umskim1>=1) .AND. (vmskjp1+vmskjm1==0)) THEN 
     1385                     risfdep(ji,jj)=0.0 ; misfdep(ji,jj) = 1 
     1386                     bathy  (ji,jj)=0.0 ; mbathy (ji,jj) = 0 
     1387                  END IF 
     1388               ! meridionnal channel and single bay 
     1389                  IF ((vmskjp1+vmskjm1>=1) .AND. (umskip1+umskim1==0)) THEN 
     1390                     risfdep(ji,jj)=0.0 ; misfdep(ji,jj) = 1 
     1391                     bathy  (ji,jj)=0.0 ; mbathy (ji,jj) = 0 
     1392                  END IF 
     1393               END IF 
     1394            END DO 
     1395         END DO 
     1396         ! ensure halo correct  
     1397         IF( lk_mpp ) THEN 
     1398            zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) ) 
     1399            zdummy(:,:) = FLOAT( mbathy(:,:)  ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:)  = INT( zdummy(:,:) ) 
     1400            CALL lbc_lnk( risfdep, 'T', 1. ) 
     1401            CALL lbc_lnk( bathy  , 'T', 1. ) 
     1402         ENDIF        
     1403      END IF 
     1404      ! 
     1405      ! 1.4 Closed single isolated point 
     1406      zmask(:,:)=1._wp 
     1407      WHERE (mbathy(:,:)==0.0) 
     1408         zmask(:,:)=0._wp 
     1409      END WHERE 
     1410      DO jj = 2, jpjm1 
     1411         DO ji = 2, jpim1 
     1412            IF(  (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN 
     1413               umskip1=zmask(ji,jj)*zmask(ji+1,jj  )  ! 1 = connexion 
     1414               vmskjp1=zmask(ji,jj)*zmask(ji  ,jj+1)  ! 1 = connexion 
     1415               umskim1=zmask(ji,jj)*zmask(ji-1,jj  )  ! 1 = connexion 
     1416               vmskjm1=zmask(ji,jj)*zmask(ji  ,jj-1)  ! 1 = connexion 
     1417               ! remove single point 
     1418               IF (umskip1+umskim1+vmskjp1+vmskjm1==0.0_wp) THEN 
     1419                  risfdep(ji,jj)=0.0 ; misfdep(ji,jj) = 1 
     1420                  bathy  (ji,jj)=0.0 ; mbathy (ji,jj) = 0 
     1421               END IF 
     1422            END IF 
     1423         END DO 
     1424      END DO 
     1425      ! 
     1426      ! ensure halo correct  
     1427      IF( lk_mpp ) THEN 
     1428         zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) ) 
     1429         zdummy(:,:) = FLOAT( mbathy(:,:)  ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:)  = INT( zdummy(:,:) ) 
     1430         CALL lbc_lnk( risfdep, 'T', 1. ) 
     1431         CALL lbc_lnk( bathy  , 'T', 1. ) 
     1432      ENDIF 
     1433      ! 
     1434      ! 2 Compute misfdep for ocean points (i.e. first wet level)  
    12981435      ! find the first ocean level such that the first level thickness  
    12991436      ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where  
    13001437      ! e3t_0 is the reference level thickness  
     1438      ! 
     1439      ! 2.0 set misfdep to 1 on open water and initialisation beneath ice shelf 
     1440      WHERE( risfdep(:,:) == 0.0 ) ;   misfdep(:,:) = 1   ! open water or grounded : set misfdep to 1   
     1441      ELSEWHERE                    ;   misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level  
     1442      END WHERE   
     1443      ! 
    13011444      DO jk = 2, jpkm1  
    13021445         zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )  
    1303          WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth )   misfdep(:,:) = jk+1  
     1446         WHERE( risfdep(:,:) > 0.0 .AND. risfdep(:,:) >= zdepth )   misfdep(:,:) = jk+1  
    13041447      END DO  
    1305       WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp) 
    1306          risfdep(:,:) = 0. ; misfdep(:,:) = 1 
    1307       END WHERE 
    1308   
    1309 ! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation 
    1310       icompt = 0  
    1311 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 
    1312       DO jl = 1, 10      
    1313          WHERE (bathy(:,:) .EQ. risfdep(:,:) ) 
    1314             misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 
    1315             mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp 
     1448      ! 
     1449      ! 2.1 fill isolated grid point in the bathymetry 
     1450      ! need to be done here to fix compatibility (will be done again in zgr_bat_ctl) 
     1451      DO jj = 2, jpjm1 
     1452         DO ji = 2, jpim1 
     1453            ibtest = MAX(  mbathy(ji-1,jj), mbathy(ji+1,jj),   & 
     1454               &           mbathy(ji,jj-1), mbathy(ji,jj+1)  ) 
     1455               IF( ibtest < mbathy(ji,jj) ) THEN 
     1456                  mbathy(ji,jj) = ibtest 
     1457                  icompt = icompt + 1 
     1458               END IF 
     1459         END DO 
     1460      END DO 
     1461      ! 
     1462      ! ensure halo correct  
     1463      zdummy(:,:) = FLOAT( mbathy(:,:) ) ; CALL lbc_lnk( zdummy, 'T', 1._wp ) ; mbathy(:,:) = INT( zdummy(:,:) ) 
     1464      ! 
     1465      IF( lk_mpp )   CALL mpp_sum( icompt ) 
     1466      IF( icompt == 0 ) THEN 
     1467         IF(lwp) WRITE(numout,*)'     no isolated ocean grid points' 
     1468      ELSE 
     1469         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points suppressed' 
     1470      ENDIF 
     1471      ! 
     1472      ! 2.2 be sure misfdep not below mbathy  
     1473      ! warning because of condition 2.0 and 2.1 we could have 'wet cell with misfdep below mbathy 
     1474      ! risfdep of these cells will be fix later on 
     1475      WHERE( misfdep > mbathy ) misfdep(:,:) = MAX( 1, mbathy(:,:) ) 
     1476      ! 
     1477      ! 3.0 Assure 2 wet cells in the column at T point and along the edge. 
     1478      ! first do T, then loop 4 times on U-1, U, V-1, V and then at T point 
     1479      ! If un-luky, digging cell U-1 can trigger case for U+1, then V-1, then V+1 
     1480      ! 
     1481      ! find the deepest isfdep level that fit the 2 wet cell on the water column 
     1482      ! on all the sides (still need 4 pass) 
     1483      DO jn = 1, 4 
     1484         zmisfdep = misfdep 
     1485         DO jj = 2, jpjm1 
     1486            DO ji = 2, jpim1 
     1487               ! ISF cell only 
     1488               IF(  (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN 
     1489                  zmbathyij  = mbathy(ji  ,jj) 
     1490                  ! set ground edge value to jpk to skip it later on 
     1491                  zmbathyip1 = mbathy(ji+1,jj) ; IF (zmbathyip1 < misfdep(ji,jj)) zmbathyip1 = jpk ! not wet cell in common on this edge 
     1492                  zmbathyim1 = mbathy(ji-1,jj) ; IF (zmbathyim1 < misfdep(ji,jj)) zmbathyim1 = jpk ! not wet cell in common on this edge 
     1493                  zmbathyjp1 = mbathy(ji,jj+1) ; IF (zmbathyjp1 < misfdep(ji,jj)) zmbathyjp1 = jpk ! not wet cell in common on this edge 
     1494                  zmbathyjm1 = mbathy(ji,jj-1) ; IF (zmbathyjm1 < misfdep(ji,jj)) zmbathyjm1 = jpk ! not wet cell in common on this edge 
     1495                  ! shallowest bathy level 
     1496                  zmbathyij = MIN(zmbathyij,zmbathyip1,zmbathyim1,zmbathyjp1,zmbathyjm1) 
     1497                  ! misfdep need to be <= zmbathyij-1 to fit 2 wet cell on the 
     1498                  ! water column 
     1499                  jk = MIN(misfdep(ji,jj),zmbathyij-1) 
     1500                  ! update isf draft if needed 
     1501                  IF ( jk < misfdep(ji,jj) ) THEN 
     1502                     zmisfdep(ji,jj) = jk 
     1503                     risfdep(ji,jj)  = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
     1504                  ELSE 
     1505                     ! isf depth not modify, nothing to do 
     1506                  END IF 
     1507               ENDIF 
     1508            END DO 
     1509         END DO 
     1510         misfdep=zmisfdep 
     1511         ! ensure halo correct before new pass  
     1512         IF( lk_mpp ) THEN 
     1513            zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) ) 
     1514            zdummy(:,:) = FLOAT( mbathy(:,:)  ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:)  = INT( zdummy(:,:) ) 
     1515            CALL lbc_lnk( risfdep, 'T', 1. ) 
     1516            CALL lbc_lnk( bathy  , 'T', 1. ) 
     1517         ENDIF 
     1518      END DO ! jn 
     1519      ! 
     1520      ! 3.1 condition block to assure connectivity every where beneath an ice shelf 
     1521      IF (ln_isfconnect) THEN 
     1522         zmask(:,:)=1.0 
     1523         WHERE (mbathy==0) 
     1524            zmask(:,:)=jpk 
    13161525         END WHERE 
    1317          WHERE (mbathy(:,:) <= 0)  
    1318             misfdep(:,:) = 0; risfdep(:,:) = 0._wp  
    1319             mbathy (:,:) = 0; bathy  (:,:) = 0._wp 
    1320          ENDWHERE 
    1321          IF( lk_mpp ) THEN 
    1322             zbathy(:,:) = FLOAT( misfdep(:,:) ) 
    1323             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1324             misfdep(:,:) = INT( zbathy(:,:) ) 
    1325             CALL lbc_lnk( risfdep, 'T', 1. ) 
    1326             CALL lbc_lnk( bathy, 'T', 1. ) 
    1327             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1328             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1329             mbathy(:,:) = INT( zbathy(:,:) ) 
    1330          ENDIF 
    1331          IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN  
    1332             misfdep( 1 ,:) = misfdep(jpim1,:)           ! local domain is cyclic east-west  
    1333             misfdep(jpi,:) = misfdep(  2  ,:)  
    1334          ENDIF 
    1335  
    1336          IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN 
    1337             mbathy( 1 ,:) = mbathy(jpim1,:)             ! local domain is cyclic east-west 
    1338             mbathy(jpi,:) = mbathy(  2  ,:) 
    1339          ENDIF 
    1340  
    1341          ! split last cell if possible (only where water column is 2 cell or less) 
    1342          DO jk = jpkm1, 1, -1 
    1343             zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
    1344             WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 
    1345                mbathy(:,:) = jk 
    1346                bathy(:,:)  = zmax 
    1347             END WHERE 
    1348          END DO 
    1349   
    1350          ! split top cell if possible (only where water column is 2 cell or less) 
    1351          DO jk = 2, jpkm1 
    1352             zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
    1353             WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy) 
    1354                misfdep(:,:) = jk 
    1355                risfdep(:,:) = zmax 
    1356             END WHERE 
    1357          END DO 
    1358  
    1359   
    1360  ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 
    1361          DO jj = 1, jpj 
    1362             DO ji = 1, jpi 
    1363                ! find the minimum change option: 
    1364                ! test bathy 
    1365                IF (risfdep(ji,jj) .GT. 1) THEN 
    1366                zbathydiff =ABS(bathy(ji,jj)   - (gdepw_1d(mbathy (ji,jj)+1) & 
    1367                  &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 
    1368                zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) & 
    1369                  &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
    1370   
    1371                   IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN 
    1372                      IF (zbathydiff .LE. zrisfdepdiff) THEN 
    1373                         bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 
    1374                         mbathy(ji,jj)= mbathy(ji,jj) + 1 
    1375                      ELSE 
    1376                         risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 
    1377                         misfdep(ji,jj) = misfdep(ji,jj) - 1 
    1378                      END IF 
     1526 
     1527         zmbathy  = mbathy 
     1528         zmisfdep = misfdep 
     1529         zrisfdep = risfdep 
     1530         WHERE (mbathy == 0) zmbathy=jpk 
     1531         DO jj = 2, jpjm1 
     1532            DO ji = 2, jpim1 
     1533               IF(  (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN 
     1534! what it should be 
     1535                  umskip1=zmask(ji,jj)*zmask(ji+1,jj  )  ! 1 = should be connected 
     1536                  umskim1=zmask(ji,jj)*zmask(ji-1,jj  )  ! 1 = should be connected 
     1537                  vmskjp1=zmask(ji,jj)*zmask(ji  ,jj+1)  ! 1 = should be connected 
     1538                  vmskjm1=zmask(ji,jj)*zmask(ji  ,jj-1)  ! 1 = should be connected 
     1539! what it is 
     1540                  umskip1_r=jpk ; umskim1_r=jpk; vmskjp1_r=jpk; vmskjm1_r=jpk 
     1541                  IF (misfdep(ji,jj) > zmbathy(ji+1,jj  )) umskip1_r=1.0 ! 1 = no effective connection 
     1542                  IF (misfdep(ji,jj) > zmbathy(ji-1,jj  )) umskim1_r=1.0 
     1543                  IF (misfdep(ji,jj) > zmbathy(ji  ,jj+1)) vmskjp1_r=1.0 
     1544                  IF (misfdep(ji,jj) > zmbathy(ji  ,jj-1)) vmskjm1_r=1.0 
     1545                  ! defining level need for connectivity 
     1546                  jk=NINT(MIN(zmbathy(ji+1,jj  ) * umskip1_r * umskip1, & 
     1547                     &        zmbathy(ji-1,jj  ) * umskim1_r * umskim1, & 
     1548                     &        zmbathy(ji  ,jj+1) * vmskjp1_r * vmskjp1, & 
     1549                     &        zmbathy(ji  ,jj-1) * vmskjm1_r * vmskjm1, & 
     1550                     &        jpk+1.0 )) ! add jpk in the MIN to avoid out of boundary later on 
     1551! if connectivity is OK or no connection needed (grounding line) or grounded, zmisfdep=misfdep 
     1552                  zmisfdep(ji,jj)=MIN(misfdep(ji,jj),jk-1) 
     1553                  ! 
     1554                  ! update isf draft if needed (need to be done here because next test is in meter and level) 
     1555                  IF ( zmisfdep(ji,jj) < misfdep(ji,jj) ) THEN 
     1556                     jk = zmisfdep(ji,jj) 
     1557                     zrisfdep(ji,jj)  = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
     1558                  END IF 
     1559    
     1560                  ! sanity check 
     1561                  ! check if we dig more than nn_kisfmax level or reach the surface 
     1562                  ! check if we dig more than rn_zisfmax meter 
     1563                  ! => if this is the case, undo what has been done before 
     1564                  IF (      (misfdep(ji,jj)-zmisfdep(ji,jj) > MIN(nn_kisfmax,misfdep(ji,jj)-2)) & 
     1565                     & .OR. (risfdep(ji,jj)-zrisfdep(ji,jj) > MIN(rn_zisfmax,risfdep(ji,jj)  )) ) THEN 
     1566                     zmisfdep(ji,jj)=misfdep(ji,jj)  
     1567                     zrisfdep(ji,jj)=risfdep(ji,jj) 
    13791568                  END IF 
    13801569               END IF 
    13811570            END DO 
    13821571         END DO 
    1383   
    1384           ! At least 2 levels for water thickness at T, U, and V point. 
    1385          DO jj = 1, jpj 
    1386             DO ji = 1, jpi 
    1387                ! find the minimum change option: 
    1388                ! test bathy 
    1389                IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
    1390                   zbathydiff  =ABS(bathy(ji,jj)  - (gdepw_1d(mbathy (ji,jj)+1)& 
    1391                     & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 
    1392                   zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  )  & 
    1393                     & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
    1394                   IF (zbathydiff .LE. zrisfdepdiff) THEN 
    1395                      mbathy(ji,jj) = mbathy(ji,jj) + 1 
    1396                      bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 
    1397                   ELSE 
    1398                      misfdep(ji,jj)= misfdep(ji,jj) - 1 
    1399                      risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 
    1400                   END IF 
    1401                ENDIF 
    1402             END DO 
    1403          END DO 
    1404   
    1405  ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1)  
    1406          DO jj = 1, jpjm1 
    1407             DO ji = 1, jpim1 
    1408                IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
    1409                   zbathydiff =ABS(bathy(ji,jj    ) - (gdepw_1d(mbathy (ji,jj)+1)   & 
    1410                     &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj  )+1)*e3zps_rat ))) 
    1411                   zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1))  & 
    1412                     &  - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 
    1413                   IF (zbathydiff .LE. zrisfdepdiff) THEN 
    1414                      mbathy(ji,jj) = mbathy(ji,jj) + 1 
    1415                      bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj  )) & 
    1416                    &    + MIN( e3zps_min, e3t_1d(mbathy(ji,jj   )+1)*e3zps_rat ) 
    1417                   ELSE 
    1418                      misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1 
    1419                      risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) & 
    1420                    &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 
    1421                   END IF 
    1422                ENDIF 
    1423             END DO 
    1424          END DO 
    1425   
     1572         misfdep=zmisfdep 
     1573         risfdep=zrisfdep 
     1574         ! ensure halo correct  
    14261575         IF( lk_mpp ) THEN 
    1427             zbathy(:,:) = FLOAT( misfdep(:,:) ) 
    1428             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1429             misfdep(:,:) = INT( zbathy(:,:) ) 
     1576            zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) ) 
     1577            zdummy(:,:) = FLOAT( mbathy(:,:)  ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:)  = INT( zdummy(:,:) ) 
    14301578            CALL lbc_lnk( risfdep, 'T', 1. ) 
    1431             CALL lbc_lnk( bathy, 'T', 1. ) 
    1432             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1433             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1434             mbathy(:,:) = INT( zbathy(:,:) ) 
     1579            CALL lbc_lnk( bathy  , 'T', 1. ) 
    14351580         ENDIF 
    1436  ! point V misdep(ji,jj) EQ mbathy(ji,jj+1)  
    1437          DO jj = 1, jpjm1 
    1438             DO ji = 1, jpim1 
    1439                IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN 
    1440                   zbathydiff =ABS(  bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) & 
    1441                    &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 
    1442                   zrisfdepdiff=ABS(risfdep(ji,jj  ) - (gdepw_1d(misfdep(ji,jj  )  )  & 
    1443                    &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )-1)*e3zps_rat ))) 
    1444                   IF (zbathydiff .LE. zrisfdepdiff) THEN 
    1445                      mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 
    1446                      bathy  (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 
    1447                   ELSE 
    1448                      misfdep(ji,jj)   = misfdep(ji,jj) - 1 
    1449                      risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat ) 
    1450                   END IF 
    1451                ENDIF 
    1452             END DO 
    1453          END DO 
    1454   
    1455   
    1456          IF( lk_mpp ) THEN  
    1457             zbathy(:,:) = FLOAT( misfdep(:,:) )  
    1458             CALL lbc_lnk( zbathy, 'T', 1. )  
    1459             misfdep(:,:) = INT( zbathy(:,:) )  
    1460             CALL lbc_lnk( risfdep, 'T', 1. )  
    1461             CALL lbc_lnk( bathy, 'T', 1. ) 
    1462             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1463             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1464             mbathy(:,:) = INT( zbathy(:,:) ) 
    1465          ENDIF  
    1466   
    1467  ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1)  
    1468          DO jj = 1, jpjm1 
    1469             DO ji = 1, jpim1 
    1470                IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
    1471                   zbathydiff =ABS(  bathy(ji  ,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 
    1472                     &   + MIN( e3zps_min, e3t_1d(mbathy (ji  ,jj)+1)*e3zps_rat ))) 
    1473                   zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) & 
    1474                     &  - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 
    1475                   IF (zbathydiff .LE. zrisfdepdiff) THEN 
    1476                      mbathy(ji,jj) = mbathy(ji,jj) + 1 
    1477                      bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 
    1478                   ELSE 
    1479                      misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 
    1480                      risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 
    1481                   END IF 
    1482                ENDIF 
     1581      END IF 
     1582      ! 
     1583      IF (ln_isfcheminey) THEN 
     1584      ! 3.2 fill hole in ice shelf (ie cell with no velocity point) 
     1585      !      => misfdep = MIN(misfdep at U, U-1, V, V-1) 
     1586      !         risfdep = gdepw(misfdep) (ie full cell) 
     1587         zmisfdep = misfdep 
     1588         WHERE (mbathy == 0) zmisfdep=jpk   ! grounded 
     1589         DO jj = 2, jpjm1 
     1590            DO ji = 2, jpim1 
     1591               ibtest    = zmisfdep(ji  ,jj  ) 
     1592               ibtestim1 = zmisfdep(ji-1,jj  ) ; ibtestip1 = zmisfdep(ji+1,jj  ) 
     1593               ibtestjm1 = zmisfdep(ji  ,jj-1) ; ibtestjp1 = zmisfdep(ji  ,jj+1) 
     1594               ! case unconected wet cell (T point) where isfdraft(ii,jj) deeper 
     1595               ! than bathy U-1/U+1/V-1/V+1 
     1596               IF( zmisfdep(ji,jj) > mbathy(ji-1,jj  ) ) ibtestim1 = jpk ! mask at U-1 
     1597               IF( zmisfdep(ji,jj) > mbathy(ji+1,jj  ) ) ibtestip1 = jpk ! mask at U+1 
     1598               IF( zmisfdep(ji,jj) > mbathy(ji  ,jj-1) ) ibtestjm1 = jpk ! mask at V-1 
     1599               IF( zmisfdep(ji,jj) > mbathy(ji  ,jj+1) ) ibtestjp1 = jpk ! mask at V+1 
     1600               ! case unconected wet cell (T point) where bathy(ii,jj) shallower 
     1601               ! than isfdraft U-1/U+1/V-1/V+1 
     1602               IF( mbathy(ji,jj) < zmisfdep(ji-1,jj  ) ) ibtestim1 = jpk ! mask at U-1 
     1603               IF( mbathy(ji,jj) < zmisfdep(ji+1,jj  ) ) ibtestip1 = jpk ! mask at U+1 
     1604               IF( mbathy(ji,jj) < zmisfdep(ji  ,jj-1) ) ibtestjm1 = jpk ! mask at V-1 
     1605               IF( mbathy(ji,jj) < zmisfdep(ji  ,jj+1) ) ibtestjp1 = jpk ! mask at V+1 
     1606               ! if surround in fully mask => mask cell  (1) 
     1607               ! if hole in the ice shelf, set to the min of surrounding (the MIN is doing the job)  
     1608               ! if misfdep is not changed, nothing is done (the MAX is doing the 
     1609               ! job) 
     1610               ibtest = MAX(ibtest,MIN(ibtestim1, ibtestip1, ibtestjm1,ibtestjp1)) 
     1611               IF (misfdep(ji,jj) < ibtest) THEN 
     1612                  misfdep(ji,jj) = ibtest 
     1613                  risfdep(ji,jj) = gdepw_1d(ibtest) 
     1614               END IF 
    14831615            ENDDO 
    14841616         ENDDO 
    1485   
    1486          IF( lk_mpp ) THEN  
    1487             zbathy(:,:) = FLOAT( misfdep(:,:) )  
    1488             CALL lbc_lnk( zbathy, 'T', 1. )  
    1489             misfdep(:,:) = INT( zbathy(:,:) )  
    1490             CALL lbc_lnk( risfdep, 'T', 1. )  
    1491             CALL lbc_lnk( bathy, 'T', 1. ) 
    1492             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1493             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1494             mbathy(:,:) = INT( zbathy(:,:) ) 
    1495          ENDIF  
    1496   
    1497  ! point U misfdep(ji,jj) EQ bathy(ji,jj+1)  
    1498          DO jj = 1, jpjm1 
    1499             DO ji = 1, jpim1 
    1500                IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
    1501                   zbathydiff =ABS(  bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) & 
    1502                       &   + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 
    1503                   zrisfdepdiff=ABS(risfdep(ji  ,jj) - (gdepw_1d(misfdep(ji  ,jj)  )  & 
    1504                       &  - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)-1)*e3zps_rat ))) 
    1505                   IF (zbathydiff .LE. zrisfdepdiff) THEN 
    1506                      mbathy(ji+1,jj)  = mbathy (ji+1,jj) + 1 
    1507                      bathy (ji+1,jj)  = gdepw_1d(mbathy (ji+1,jj)  )  & 
    1508                       &   + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 
    1509                   ELSE 
    1510                      misfdep(ji,jj)   = misfdep(ji  ,jj) - 1 
    1511                      risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) & 
    1512                       &   - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat ) 
    1513                   END IF 
    1514                ENDIF 
    1515             ENDDO 
    1516          ENDDO 
    1517   
     1617         WHERE( misfdep == jpk)   ! ground case (1) 
     1618            mbathy = 0 ; bathy = 0.0_wp ; misfdep = 1 ; risfdep = 0.0_wp 
     1619         END WHERE 
     1620         ! 
     1621         ! ensure halo correct 
    15181622         IF( lk_mpp ) THEN 
    1519             zbathy(:,:) = FLOAT( misfdep(:,:) ) 
    1520             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1521             misfdep(:,:) = INT( zbathy(:,:) ) 
     1623            zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) ) 
     1624            zdummy(:,:) = FLOAT( mbathy(:,:)  ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:)  = INT( zdummy(:,:) ) 
    15221625            CALL lbc_lnk( risfdep, 'T', 1. ) 
    1523             CALL lbc_lnk( bathy, 'T', 1. ) 
    1524             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1525             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1526             mbathy(:,:) = INT( zbathy(:,:) ) 
    1527          ENDIF 
    1528       END DO 
    1529       ! end dig bathy/ice shelf to be compatible 
    1530       ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness 
    1531       DO jl = 1,20 
    1532   
    1533  ! remove single point "bay" on isf coast line in the ice shelf draft' 
    1534          DO jk = 2, jpk 
    1535             WHERE (misfdep==0) misfdep=jpk 
    1536             zmask=0 
    1537             WHERE (misfdep .LE. jk) zmask=1 
    1538             DO jj = 2, jpjm1 
    1539                DO ji = 2, jpim1 
    1540                   IF (misfdep(ji,jj) .EQ. jk) THEN 
    1541                      ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 
    1542                      IF (ibtest .LE. 1) THEN 
    1543                         risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1 
    1544                         IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) misfdep(ji,jj) = jpk 
    1545                      END IF 
    1546                   END IF 
    1547                END DO 
    1548             END DO 
    1549          END DO 
    1550          WHERE (misfdep==jpk) 
    1551              misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0. 
    1552          END WHERE 
    1553          IF( lk_mpp ) THEN 
    1554             zbathy(:,:) = FLOAT( misfdep(:,:) ) 
    1555             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1556             misfdep(:,:) = INT( zbathy(:,:) ) 
    1557             CALL lbc_lnk( risfdep, 'T', 1. ) 
    1558             CALL lbc_lnk( bathy, 'T', 1. ) 
    1559             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1560             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1561             mbathy(:,:) = INT( zbathy(:,:) ) 
    1562          ENDIF 
    1563   
    1564  ! remove single point "bay" on bathy coast line beneath an ice shelf' 
    1565          DO jk = jpk,1,-1 
    1566             zmask=0 
    1567             WHERE (mbathy .GE. jk ) zmask=1 
    1568             DO jj = 2, jpjm1 
    1569                DO ji = 2, jpim1 
    1570                   IF (mbathy(ji,jj) .EQ. jk .AND. misfdep(ji,jj) .GE. 2) THEN 
    1571                      ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 
    1572                      IF (ibtest .LE. 1) THEN 
    1573                         bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1 
    1574                         IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) mbathy(ji,jj) = 0 
    1575                      END IF 
    1576                   END IF 
    1577                END DO 
    1578             END DO 
    1579          END DO 
    1580          WHERE (mbathy==0) 
    1581              misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0. 
    1582          END WHERE 
    1583          IF( lk_mpp ) THEN 
    1584             zbathy(:,:) = FLOAT( misfdep(:,:) ) 
    1585             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1586             misfdep(:,:) = INT( zbathy(:,:) ) 
    1587             CALL lbc_lnk( risfdep, 'T', 1. ) 
    1588             CALL lbc_lnk( bathy, 'T', 1. ) 
    1589             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1590             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1591             mbathy(:,:) = INT( zbathy(:,:) ) 
    1592          ENDIF 
    1593   
    1594  ! fill hole in ice shelf 
    1595          zmisfdep = misfdep 
    1596          zrisfdep = risfdep 
    1597          WHERE (zmisfdep .LE. 1) zmisfdep=jpk 
    1598          DO jj = 2, jpjm1 
    1599             DO ji = 2, jpim1 
    1600                ibtestim1 = zmisfdep(ji-1,jj  ) ; ibtestip1 = zmisfdep(ji+1,jj  ) 
    1601                ibtestjm1 = zmisfdep(ji  ,jj-1) ; ibtestjp1 = zmisfdep(ji  ,jj+1) 
    1602                IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj  ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj  ) - 1) 
    1603                IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj  ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj  ) - 1) 
    1604                IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji  ,jj-1) - 1) 
    1605                IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji  ,jj+1) - 1) 
    1606                ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 
    1607                IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN 
    1608                   mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp 
    1609                END IF 
    1610                IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) .GE. 2) THEN 
    1611                   misfdep(ji,jj) = ibtest 
    1612                   risfdep(ji,jj) = gdepw_1d(ibtest) 
    1613                ENDIF 
    1614             ENDDO 
    1615          ENDDO 
    1616   
    1617          IF( lk_mpp ) THEN  
    1618             zbathy(:,:) = FLOAT( misfdep(:,:) )  
    1619             CALL lbc_lnk( zbathy, 'T', 1. )  
    1620             misfdep(:,:) = INT( zbathy(:,:) )  
    1621             CALL lbc_lnk( risfdep, 'T', 1. )  
    1622             CALL lbc_lnk( bathy, 'T', 1. ) 
    1623             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1624             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1625             mbathy(:,:) = INT( zbathy(:,:) ) 
    1626          ENDIF  
    1627  ! 
    1628  !! fill hole in bathymetry 
    1629          zmbathy (:,:)=mbathy (:,:) 
    1630          DO jj = 2, jpjm1 
    1631             DO ji = 2, jpim1 
    1632                ibtestim1 = zmbathy(ji-1,jj  ) ; ibtestip1 = zmbathy(ji+1,jj  ) 
    1633                ibtestjm1 = zmbathy(ji  ,jj-1) ; ibtestjp1 = zmbathy(ji  ,jj+1) 
    1634                IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj  ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj  ) + 1) 
    1635                IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj  ) ) ibtestip1 = 0 
    1636                IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj-1) ) ibtestjm1 = 0 
    1637                IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj+1) ) ibtestjp1 = 0 
    1638                ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 
    1639                IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN 
    1640                   mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 
    1641                END IF 
    1642                IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) .GE. 2) THEN 
    1643                   mbathy(ji,jj) = ibtest 
    1644                   bathy(ji,jj)  = gdepw_1d(ibtest+1)  
    1645                ENDIF 
    1646             END DO 
    1647          END DO 
    1648          IF( lk_mpp ) THEN  
    1649             zbathy(:,:) = FLOAT( misfdep(:,:) )  
    1650             CALL lbc_lnk( zbathy, 'T', 1. )  
    1651             misfdep(:,:) = INT( zbathy(:,:) )  
    1652             CALL lbc_lnk( risfdep, 'T', 1. )  
    1653             CALL lbc_lnk( bathy, 'T', 1. ) 
    1654             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1655             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1656             mbathy(:,:) = INT( zbathy(:,:) ) 
    1657          ENDIF  
    1658  ! if not compatible after all check (ie U point water column less than 2 cells), mask U 
    1659          DO jj = 1, jpjm1 
    1660             DO ji = 1, jpim1 
    1661                IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN 
    1662                   mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ; 
    1663                END IF 
    1664             END DO 
    1665          END DO 
    1666          IF( lk_mpp ) THEN  
    1667             zbathy(:,:) = FLOAT( misfdep(:,:) )  
    1668             CALL lbc_lnk( zbathy, 'T', 1. )  
    1669             misfdep(:,:) = INT( zbathy(:,:) )  
    1670             CALL lbc_lnk( risfdep, 'T', 1. )  
    1671             CALL lbc_lnk( bathy, 'T', 1. ) 
    1672             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1673             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1674             mbathy(:,:) = INT( zbathy(:,:) ) 
    1675          ENDIF  
    1676  ! if not compatible after all check (ie U point water column less than 2 cells), mask U 
    1677          DO jj = 1, jpjm1 
    1678             DO ji = 1, jpim1 
    1679                IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN 
    1680                   mbathy(ji+1,jj)  = mbathy(ji+1,jj) - 1;   bathy(ji+1,jj)   = gdepw_1d(mbathy(ji+1,jj)+1) ; 
    1681                END IF 
    1682             END DO 
    1683          END DO 
    1684          IF( lk_mpp ) THEN  
    1685             zbathy(:,:) = FLOAT( misfdep(:,:) )  
    1686             CALL lbc_lnk( zbathy, 'T', 1. )  
    1687             misfdep(:,:) = INT( zbathy(:,:) )  
    1688             CALL lbc_lnk( risfdep, 'T', 1. )  
    1689             CALL lbc_lnk( bathy, 'T', 1. ) 
    1690             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1691             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1692             mbathy(:,:) = INT( zbathy(:,:) ) 
    1693          ENDIF  
    1694  ! if not compatible after all check (ie V point water column less than 2 cells), mask V 
    1695          DO jj = 1, jpjm1 
    1696             DO ji = 1, jpi 
    1697                IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN 
    1698                   mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ; 
    1699                END IF 
    1700             END DO 
    1701          END DO 
    1702          IF( lk_mpp ) THEN  
    1703             zbathy(:,:) = FLOAT( misfdep(:,:) )  
    1704             CALL lbc_lnk( zbathy, 'T', 1. )  
    1705             misfdep(:,:) = INT( zbathy(:,:) )  
    1706             CALL lbc_lnk( risfdep, 'T', 1. )  
    1707             CALL lbc_lnk( bathy, 'T', 1. ) 
    1708             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1709             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1710             mbathy(:,:) = INT( zbathy(:,:) ) 
    1711          ENDIF  
    1712  ! if not compatible after all check (ie V point water column less than 2 cells), mask V 
    1713          DO jj = 1, jpjm1 
    1714             DO ji = 1, jpi 
    1715                IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN 
    1716                   mbathy(ji,jj+1)  = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1)   = gdepw_1d(mbathy(ji,jj+1)+1) ; 
    1717                END IF 
    1718             END DO 
    1719          END DO 
    1720          IF( lk_mpp ) THEN  
    1721             zbathy(:,:) = FLOAT( misfdep(:,:) )  
    1722             CALL lbc_lnk( zbathy, 'T', 1. )  
    1723             misfdep(:,:) = INT( zbathy(:,:) )  
    1724             CALL lbc_lnk( risfdep, 'T', 1. )  
    1725             CALL lbc_lnk( bathy, 'T', 1. ) 
    1726             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1727             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1728             mbathy(:,:) = INT( zbathy(:,:) ) 
    1729          ENDIF  
    1730  ! if not compatible after all check, mask T 
    1731          DO jj = 1, jpj 
    1732             DO ji = 1, jpi 
    1733                IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN 
    1734                   misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj)  = 0 ; bathy(ji,jj)   = 0._wp ; 
    1735                END IF 
    1736             END DO 
    1737          END DO 
    1738   
    1739          WHERE (mbathy(:,:) == 1) 
    1740             mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp 
    1741          END WHERE 
    1742       END DO  
    1743 ! end check compatibility ice shelf/bathy 
    1744       ! remove very shallow ice shelf (less than ~ 10m if 75L) 
    1745       WHERE (misfdep(:,:) <= 5) 
    1746          misfdep = 1; risfdep = 0.0_wp; 
    1747       END WHERE 
    1748  
    1749       IF( icompt == 0 ) THEN  
    1750          IF(lwp) WRITE(numout,*)'     no points with ice shelf too close to bathymetry'  
    1751       ELSE  
    1752          IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry'  
    1753       ENDIF  
    1754  
    1755       CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 
    1756       CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 
    1757  
     1626            CALL lbc_lnk( bathy  , 'T', 1. ) 
     1627         END IF 
     1628      END IF 
     1629      ! 
     1630      CALL wrk_dealloc( jpi, jpj, zmask, zdummy, zrisfdep) 
     1631      CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy) 
     1632      ! 
    17581633      IF( nn_timing == 1 )  CALL timing_stop('zgr_isf') 
    17591634       
Note: See TracChangeset for help on using the changeset viewer.