Changeset 9316


Ignore:
Timestamp:
2018-02-09T13:07:18+01:00 (2 years ago)
Author:
mathiot
Message:

update domzgr and namelist (rewrite of zgr_isf based on discussion with P. Holland, R. Smith, A. Siaahan and J. Gregory)

Location:
branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/CONFIG/SHARED/namelist_ref

    r8447 r9316  
    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_isfhw_deep    = 1.e-3       ! minimum water column thickness to declare a column 'wet' 
     124   rn_isfshallow    = 0.0         ! bathy between surface and rn_isfshallow are declared shallow 
     125   rn_isfhw_shallow = 1.e-3       ! minimum water column thickness to declare a column 'wet' in case of shallow isf 
     126   ln_isfchannel    = .false.     ! remove channel (based on 2d mask build from isfdraft-bathy) 
     127   ln_isfconnect    = .false.     ! force connection under the ice shelf (based on 2d mask build from isfdraft-bathy) 
     128   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) 
     129   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) 
     130   ln_isfcheminey   = .false.     ! close cheminey 
    118131/ 
    119132!----------------------------------------------------------------------- 
  • branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r6487 r9316  
    12631263      !!                Bathymetry and isfdraft are modified (dig/close) to respect 
    12641264      !!                this criterion. 
    1265       !!                  
     1265      !!                digging and filling base on depth criterion only 
     1266      !!                          1.0 = set iceshelf to the minimum depth allowed 
     1267      !!                          1.1 = ground ice shelf if water column less than X m 
     1268      !!                          1.2 = ensure a minimum thickness for iceshelf cavity in shallow water 
     1269      !!                          1.3 = remove channels and single point 'bay' 
     1270      !!                          1.4 = close single isolated point 
     1271      !!                compute level 
     1272      !!                          2.0 = compute level 
     1273      !!                          2.1 = ensure misfdep is not below bathymetry after step 2.0 
     1274      !!                digging and filling base on level criterion only                 
     1275      !!                          3.0 = dig to fit the 2 water level criterion (closed pool possible after this step) 
     1276      !!                          3.1 = dig enough to ensure connectivity of all the cell beneath an ice shelf  
     1277      !!                                (most of the close pool should be remove after this step) 
     1278      !!                          3.2 = fill chimney 
    12661279      !!    
    12671280      !! ** Action  : - test compatibility between isfdraft and bathy  
    12681281      !!              - bathy and isfdraft are modified 
    12691282      !!---------------------------------------------------------------------- 
    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) 
     1283      !!   
     1284      INTEGER  ::   ji, jj, jk, jn 
     1285      INTEGER  ::   ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1   ! (ISF) 
     1286      INTEGER  ::   zmbathyij, zmbathyip1, zmbathyim1, zmbathyjp1, zmbathyjm1 
     1287      INTEGER , POINTER, DIMENSION(:,:)   ::   zmisfdep, zmbathy         ! 2D workspace (ISH) 
     1288      ! 
     1289      REAL(wp) ::   zdepth, vmskjp1, vmskjm1, umskip1, umskim1, umskip1_r, umskim1_r, vmskjp1_r, vmskjm1_r 
     1290      REAL(wp), POINTER, DIMENSION(:,:)   ::   zdummy, zmask, zrisfdep   ! 2D workspace (ISH) 
     1291      ! 
     1292      !!----------------- 
     1293      ! NAMELIST 
     1294      INTEGER  ::   ios 
     1295      INTEGER  :: nn_kisfmax    = 999. !   
     1296      REAL(wp) :: rn_isfdep_min = 10.0_wp ! ice shelf minimal thickness  
     1297      REAL(wp) :: rn_isfhw_deep = 1.e-3 , rn_isfhw_shallow = 1.0e-3 ! threshold to define grounding line in deep/shallow water 
     1298      REAL(wp) :: rn_isfshallow = 0.0_wp  ! threshold to define shallow ice shelf cavity 
     1299      REAL(wp) :: rn_zisfmax    = 6000.0_wp ! maximun meter of ice we are allowed to dig to assure connectivity 
     1300      LOGICAL  :: ln_isfcheminey= .FALSE. , ln_isfconnect= .FALSE. , ln_isfchannel= .FALSE. !  remove cheminey, assure connectivity 
    12841301      !!--------------------------------------------------------------------- 
     1302      NAMELIST/namzgr_isf/nn_kisfmax, rn_zisfmax, & 
     1303              &           rn_isfdep_min, rn_isfhw_deep, rn_isfhw_shallow, rn_isfshallow, & 
     1304              &           ln_isfcheminey, ln_isfconnect, ln_isfchannel 
    12851305      ! 
    12861306      IF( nn_timing == 1 )  CALL timing_start('zgr_isf') 
    12871307      ! 
    1288       CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep) 
     1308      CALL wrk_alloc( jpi, jpj, zdummy, zmask, zrisfdep) 
    12891309      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)  
     1310      ! 
     1311      ! 
     1312      REWIND( numnam_ref )              ! Namelist namzgr_isf in reference namelist : ice shelf geometry definition 
     1313      READ  ( numnam_ref, namzgr_isf, IOSTAT = ios, ERR = 901) 
     1314901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_isf in reference namelist', lwp ) 
     1315 
     1316      REWIND( numnam_cfg )              ! Namelist namzgr_sco in configuration namelist : ice shelf geometry definition 
     1317      READ  ( numnam_cfg, namzgr_isf, IOSTAT = ios, ERR = 902 ) 
     1318902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_isf in configuration namelist', lwp ) 
     1319      IF(lwm) WRITE ( numond, namzgr_isf ) 
     1320      ! 
     1321      IF (lwp) THEN 
     1322         WRITE(numout,*) ' nn_kisfmax       = ',nn_kisfmax           ! 
     1323         WRITE(numout,*) ' rn_zisfmax       = ',rn_zisfmax           ! 
     1324         WRITE(numout,*) ' rn_isfdep_min    = ',rn_isfdep_min        ! 
     1325         WRITE(numout,*) ' rn_isfhw_deep    = ',rn_isfhw_deep        ! 
     1326         WRITE(numout,*) ' rn_isfhw_shallow = ',rn_isfhw_shallow     ! 
     1327         WRITE(numout,*) ' rn_isfshallow    = ',rn_isfshallow        ! 
     1328         WRITE(numout,*) ' ln_isfcheminey   = ',ln_isfcheminey       ! 
     1329         WRITE(numout,*) ' ln_isfconnect    = ',ln_isfconnect        ! 
     1330         WRITE(numout,*) ' ln_isfchannel    = ',ln_isfchannel        ! 
     1331      END IF 
     1332      ! 
     1333      ! 1.0 set iceshelf to the minimum depth allowed 
     1334      WHERE(risfdep(:,:) > 0.0_wp .AND. risfdep(:,:) < MAX(rn_isfdep_min,e3t_1d(1))) 
     1335         risfdep(:,:)=rn_isfdep_min 
     1336      END WHERE 
     1337      !   
     1338      ! 1.1 ground ice shelf if water column less than X m 
     1339      WHERE( bathy(:,:) - risfdep(:,:) < rn_isfhw_deep )  
     1340         risfdep(:,:)=0.0 ; misfdep(:,:) = 1 
     1341         bathy  (:,:)=0.0 ; mbathy (:,:) = 0 
     1342      END WHERE 
     1343      ! 
     1344      ! 1.2 ensure a minimum thickness for iceshelf cavity in shallow water 
     1345      WHERE( (bathy(:,:) < rn_isfshallow) .AND. (bathy(:,:) - risfdep(:,:) < rn_isfhw_shallow) )  
     1346         risfdep(:,:)=0.0 ; misfdep(:,:) = 1 
     1347         bathy  (:,:)=0.0 ; mbathy (:,:) = 0 
     1348      END WHERE 
     1349      ! 
     1350      ! 1.3 Remove channels and single point 'bay'. 
     1351      ! channel could be created if connectivity is enforced later. 
     1352      IF (ln_isfchannel) THEN 
     1353         zmask(:,:)=1._wp 
     1354         WHERE (mbathy(:,:)==0.0) 
     1355            zmask(:,:)=0._wp 
     1356         END WHERE 
     1357         DO jj = 2, jpjm1 
     1358            DO ji = 2, jpim1 
     1359               IF(  (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN 
     1360                  umskip1=zmask(ji,jj)*zmask(ji+1,jj  )  ! 1 = connexion 
     1361                  vmskjp1=zmask(ji,jj)*zmask(ji  ,jj+1)  ! 1 = connexion 
     1362                  umskim1=zmask(ji,jj)*zmask(ji-1,jj  )  ! 1 = connexion 
     1363                  vmskjm1=zmask(ji,jj)*zmask(ji  ,jj-1)  ! 1 = connexion 
     1364               ! zonal channel and single bay 
     1365                  IF ((umskip1+umskim1>=1) .AND. (vmskjp1+vmskjm1==0)) THEN 
     1366                     risfdep(ji,jj)=0.0 ; misfdep(ji,jj) = 1 
     1367                     bathy  (ji,jj)=0.0 ; mbathy (ji,jj) = 0 
     1368                  END IF 
     1369               ! meridionnal channel and single bay 
     1370                  IF ((vmskjp1+vmskjm1>=1) .AND. (umskip1+umskim1==0)) THEN 
     1371                     risfdep(ji,jj)=0.0 ; misfdep(ji,jj) = 1 
     1372                     bathy  (ji,jj)=0.0 ; mbathy (ji,jj) = 0 
     1373                  END IF 
     1374               END IF 
     1375            END DO 
     1376         END DO 
     1377         ! ensure halo correct  
     1378         IF( lk_mpp ) THEN 
     1379            zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) ) 
     1380            zdummy(:,:) = FLOAT( mbathy(:,:)  ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:)  = INT( zdummy(:,:) ) 
     1381            CALL lbc_lnk( risfdep, 'T', 1. ) 
     1382            CALL lbc_lnk( bathy  , 'T', 1. ) 
     1383         ENDIF        
     1384      END IF 
     1385      ! 
     1386      ! 1.4 Closed single isolated point 
     1387      zmask(:,:)=1._wp 
     1388      WHERE (mbathy(:,:)==0.0) 
     1389         zmask(:,:)=0._wp 
     1390      END WHERE 
     1391      DO jj = 2, jpjm1 
     1392         DO ji = 2, jpim1 
     1393            IF(  (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN 
     1394               umskip1=zmask(ji,jj)*zmask(ji+1,jj  )  ! 1 = connexion 
     1395               vmskjp1=zmask(ji,jj)*zmask(ji  ,jj+1)  ! 1 = connexion 
     1396               umskim1=zmask(ji,jj)*zmask(ji-1,jj  )  ! 1 = connexion 
     1397               vmskjm1=zmask(ji,jj)*zmask(ji  ,jj-1)  ! 1 = connexion 
     1398               ! remove single point 
     1399               IF (umskip1+umskim1+vmskjp1+vmskjm1==0.0_wp) THEN 
     1400                  risfdep(ji,jj)=0.0 ; misfdep(ji,jj) = 1 
     1401                  bathy  (ji,jj)=0.0 ; mbathy (ji,jj) = 0 
     1402               END IF 
     1403            END IF 
     1404         END DO 
     1405      END DO 
     1406      ! ensure halo correct  
     1407      IF( lk_mpp ) THEN 
     1408         zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) ) 
     1409         zdummy(:,:) = FLOAT( mbathy(:,:)  ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:)  = INT( zdummy(:,:) ) 
     1410         CALL lbc_lnk( risfdep, 'T', 1. ) 
     1411         CALL lbc_lnk( bathy  , 'T', 1. ) 
     1412      ENDIF 
     1413      ! 
     1414      ! 2 Compute misfdep for ocean points (i.e. first wet level)  
    12981415      ! find the first ocean level such that the first level thickness  
    12991416      ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where  
    13001417      ! e3t_0 is the reference level thickness  
     1418      ! 
     1419      ! 2.0 set misfdep to 1 on open water and initialisation beneath ice shelf 
     1420      WHERE( risfdep(:,:) == 0.0 ) ;   misfdep(:,:) = 1   ! open water or grounded : set misfdep to 1   
     1421      ELSEWHERE                    ;   misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level  
     1422      END WHERE   
     1423      ! 
    13011424      DO jk = 2, jpkm1  
    13021425         zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )  
    1303          WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth )   misfdep(:,:) = jk+1  
     1426         WHERE( risfdep(:,:) > 0.0 .AND. risfdep(:,:) >= zdepth )   misfdep(:,:) = jk+1  
    13041427      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 
     1428      ! 
     1429      ! 2.1 be sure misfdep not below mbathy  
     1430      ! warning because of condition 4 we could have 'wet cell with misfdep below mbathy 
     1431      ! risfdep of these cells will be fix later on (see 5) 
     1432      WHERE( misfdep > mbathy .AND. mbathy > 0 ) misfdep=mbathy 
     1433      ! 
     1434      ! 3.0 Assure 2 wet cells in the column at T point and along the edge. 
     1435      ! first do T, then loop 4 times on U-1, U, V-1, V and then at T point 
     1436      ! If un-luky, digging cell U-1 can trigger case for U+1, then V-1, then V+1 
     1437      ! 
     1438      ! find the deepest isfdep level that fit the 2 wet cell on the water column 
     1439      ! on all the sides (still need 4 pass) 
     1440      DO jn = 1, 4 
     1441         zmisfdep = misfdep 
     1442         DO jj = 2, jpjm1 
     1443            DO ji = 2, jpim1 
     1444               ! ISF cell only 
     1445               IF(  (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN 
     1446                  zmbathyij  = mbathy(ji  ,jj) 
     1447                  ! set ground edge value to jpk to skip it later on 
     1448                  zmbathyip1 = mbathy(ji+1,jj) ; IF (zmbathyip1 < misfdep(ji,jj)) zmbathyip1 = jpk ! not wet cell in common on this edge 
     1449                  zmbathyim1 = mbathy(ji-1,jj) ; IF (zmbathyim1 < misfdep(ji,jj)) zmbathyim1 = jpk ! not wet cell in common on this edge 
     1450                  zmbathyjp1 = mbathy(ji,jj+1) ; IF (zmbathyjp1 < misfdep(ji,jj)) zmbathyjp1 = jpk ! not wet cell in common on this edge 
     1451                  zmbathyjm1 = mbathy(ji,jj-1) ; IF (zmbathyjm1 < misfdep(ji,jj)) zmbathyjm1 = jpk ! not wet cell in common on this edge 
     1452                  ! shallowest bathy level 
     1453                  zmbathyij = MIN(zmbathyij,zmbathyip1,zmbathyim1,zmbathyjp1,zmbathyjm1) 
     1454                  ! misfdep need to be <= zmbathyij-1 to fit 2 wet cell on the 
     1455                  ! water column 
     1456                  jk = MIN(misfdep(ji,jj),zmbathyij-1) 
     1457                  ! update isf draft if needed 
     1458                  IF ( jk < misfdep(ji,jj) ) THEN 
     1459                     zmisfdep(ji,jj) = jk 
     1460                     risfdep(ji,jj)  = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
     1461                  ELSE 
     1462                     ! isf depth not modify, nothing to do 
     1463                  END IF 
     1464               ENDIF 
     1465            END DO 
     1466         END DO 
     1467         misfdep=zmisfdep 
     1468         ! ensure halo correct before new pass  
     1469         IF( lk_mpp ) THEN 
     1470            zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) ) 
     1471            zdummy(:,:) = FLOAT( mbathy(:,:)  ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:)  = INT( zdummy(:,:) ) 
     1472            CALL lbc_lnk( risfdep, 'T', 1. ) 
     1473            CALL lbc_lnk( bathy  , 'T', 1. ) 
     1474         ENDIF 
     1475      END DO ! jn 
     1476      ! 
     1477      ! 3.1 condition block to assure connectivity every where beneath an ice shelf 
     1478      IF (ln_isfconnect) THEN 
     1479         zmask(:,:)=1.0 
     1480         WHERE (mbathy==0) 
     1481            zmask(:,:)=jpk 
    13161482         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 
     1483 
     1484         zmbathy  = mbathy 
     1485         zmisfdep = misfdep 
     1486         zrisfdep = risfdep 
     1487         WHERE (mbathy == 0) zmbathy=jpk 
     1488         DO jj = 2, jpjm1 
     1489            DO ji = 2, jpim1 
     1490               IF(  (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN 
     1491! what it should be 
     1492                  umskip1=zmask(ji,jj)*zmask(ji+1,jj  )  ! 1 = should be connected 
     1493                  umskim1=zmask(ji,jj)*zmask(ji-1,jj  )  ! 1 = should be connected 
     1494                  vmskjp1=zmask(ji,jj)*zmask(ji  ,jj+1)  ! 1 = should be connected 
     1495                  vmskjm1=zmask(ji,jj)*zmask(ji  ,jj-1)  ! 1 = should be connected 
     1496! what it is 
     1497                  umskip1_r=jpk ; umskim1_r=jpk; vmskjp1_r=jpk; vmskjm1_r=jpk 
     1498                  IF (misfdep(ji,jj) > zmbathy(ji+1,jj  )) umskip1_r=1.0 ! 1 = no effective connection 
     1499                  IF (misfdep(ji,jj) > zmbathy(ji-1,jj  )) umskim1_r=1.0 
     1500                  IF (misfdep(ji,jj) > zmbathy(ji  ,jj+1)) vmskjp1_r=1.0 
     1501                  IF (misfdep(ji,jj) > zmbathy(ji  ,jj-1)) vmskjm1_r=1.0 
     1502                  ! defining level need for connectivity 
     1503                  jk=NINT(MIN(zmbathy(ji+1,jj  ) * umskip1_r * umskip1, & 
     1504                     &        zmbathy(ji-1,jj  ) * umskim1_r * umskim1, & 
     1505                     &        zmbathy(ji  ,jj+1) * vmskjp1_r * vmskjp1, & 
     1506                     &        zmbathy(ji  ,jj-1) * vmskjm1_r * vmskjm1, & 
     1507                     &        jpk+1.0 )) ! add jpk in the MIN to avoid out of boundary later on 
     1508! if connectivity is OK or no connection needed (grounding line) or grounded, zmisfdep=misfdep 
     1509                  zmisfdep(ji,jj)=MIN(misfdep(ji,jj),jk-1) 
     1510                  ! 
     1511                  ! update isf draft if needed (need to be done here because next test is in meter and level) 
     1512                  IF ( zmisfdep(ji,jj) < misfdep(ji,jj) ) THEN 
     1513                     jk = zmisfdep(ji,jj) 
     1514                     zrisfdep(ji,jj)  = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
     1515                  END IF 
     1516    
     1517                  ! sanity check 
     1518                  ! check if we dig more than nn_kisfmax level or reach the surface 
     1519                  ! check if we dig more than rn_zisfmax meter 
     1520                  ! => if this is the case, undo what has been done before 
     1521                  IF (      (misfdep(ji,jj)-zmisfdep(ji,jj) > MIN(nn_kisfmax,misfdep(ji,jj)-2)) & 
     1522                     & .OR. (risfdep(ji,jj)-zrisfdep(ji,jj) > MIN(rn_zisfmax,risfdep(ji,jj)  )) ) THEN 
     1523                     zmisfdep(ji,jj)=misfdep(ji,jj)  
     1524                     zrisfdep(ji,jj)=risfdep(ji,jj) 
    13791525                  END IF 
    13801526               END IF 
    13811527            END DO 
    13821528         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   
     1529         misfdep=zmisfdep 
     1530         risfdep=zrisfdep 
     1531         ! ensure halo correct  
    14261532         IF( lk_mpp ) THEN 
    1427             zbathy(:,:) = FLOAT( misfdep(:,:) ) 
    1428             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1429             misfdep(:,:) = INT( zbathy(:,:) ) 
     1533            zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) ) 
     1534            zdummy(:,:) = FLOAT( mbathy(:,:)  ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:)  = INT( zdummy(:,:) ) 
    14301535            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(:,:) ) 
     1536            CALL lbc_lnk( bathy  , 'T', 1. ) 
    14351537         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 
     1538      END IF 
     1539      ! 
     1540      IF (ln_isfcheminey) THEN 
     1541      ! 3.2 fill hole in ice shelf (ie cell with no velocity point) 
     1542      !      => misfdep = MIN(misfdep at U, U-1, V, V-1) 
     1543      !         risfdep = gdepw(misfdep) (ie full cell) 
     1544         zmisfdep = misfdep 
     1545         WHERE (mbathy == 0) zmisfdep=jpk   ! grounded 
     1546         DO jj = 2, jpjm1 
     1547            DO ji = 2, jpim1 
     1548               ibtest    = zmisfdep(ji  ,jj  ) 
     1549               ibtestim1 = zmisfdep(ji-1,jj  ) ; ibtestip1 = zmisfdep(ji+1,jj  ) 
     1550               ibtestjm1 = zmisfdep(ji  ,jj-1) ; ibtestjp1 = zmisfdep(ji  ,jj+1) 
     1551               ! case unconected wet cell (T point) where isfdraft(ii,jj) deeper 
     1552               ! than bathy U-1/U+1/V-1/V+1 
     1553               IF( zmisfdep(ji,jj) > mbathy(ji-1,jj  ) ) ibtestim1 = jpk ! mask at U-1 
     1554               IF( zmisfdep(ji,jj) > mbathy(ji+1,jj  ) ) ibtestip1 = jpk ! mask at U+1 
     1555               IF( zmisfdep(ji,jj) > mbathy(ji  ,jj-1) ) ibtestjm1 = jpk ! mask at V-1 
     1556               IF( zmisfdep(ji,jj) > mbathy(ji  ,jj+1) ) ibtestjp1 = jpk ! mask at V+1 
     1557               ! case unconected wet cell (T point) where bathy(ii,jj) shallower 
     1558               ! than isfdraft U-1/U+1/V-1/V+1 
     1559               IF( mbathy(ji,jj) < zmisfdep(ji-1,jj  ) ) ibtestim1 = jpk ! mask at U-1 
     1560               IF( mbathy(ji,jj) < zmisfdep(ji+1,jj  ) ) ibtestip1 = jpk ! mask at U+1 
     1561               IF( mbathy(ji,jj) < zmisfdep(ji  ,jj-1) ) ibtestjm1 = jpk ! mask at V-1 
     1562               IF( mbathy(ji,jj) < zmisfdep(ji  ,jj+1) ) ibtestjp1 = jpk ! mask at V+1 
     1563               ! if surround in fully mask => mask cell  (1) 
     1564               ! if hole in the ice shelf, set to the min of surrounding (the MIN is doing the job)  
     1565               ! if misfdep is not changed, nothing is done (the MAX is doing the 
     1566               ! job) 
     1567               ibtest = MAX(ibtest,MIN(ibtestim1, ibtestip1, ibtestjm1,ibtestjp1)) 
     1568               IF (misfdep(ji,jj) < ibtest) THEN 
     1569                  misfdep(ji,jj) = ibtest 
     1570                  risfdep(ji,jj) = gdepw_1d(ibtest) 
     1571               END IF 
    14831572            ENDDO 
    14841573         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   
     1574         WHERE( misfdep == jpk)   ! ground case (1) 
     1575            mbathy = 0 ; bathy = 0.0_wp ; misfdep = 1 ; risfdep = 0.0_wp 
     1576         END WHERE 
     1577         ! 
     1578         ! ensure halo correct 
    15181579         IF( lk_mpp ) THEN 
    1519             zbathy(:,:) = FLOAT( misfdep(:,:) ) 
    1520             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1521             misfdep(:,:) = INT( zbathy(:,:) ) 
     1580            zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) ) 
     1581            zdummy(:,:) = FLOAT( mbathy(:,:)  ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:)  = INT( zdummy(:,:) ) 
    15221582            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  
     1583            CALL lbc_lnk( bathy  , 'T', 1. ) 
     1584         END IF 
     1585      END IF 
     1586      ! 
     1587      CALL wrk_dealloc( jpi, jpj, zmask, zdummy, zrisfdep) 
     1588      CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy) 
     1589      ! 
    17581590      IF( nn_timing == 1 )  CALL timing_stop('zgr_isf') 
    17591591       
Note: See TracChangeset for help on using the changeset viewer.