Changeset 14573 for branches/NERC
- Timestamp:
- 2021-03-03T14:39:55+01:00 (3 years ago)
- 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 116 116 !!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above] 117 117 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 118 133 / 119 134 !----------------------------------------------------------------------- -
branches/NERC/dev_r5518_GO6_r9321_isf_merge_UKESM_AIS/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
r14572 r14573 25 25 USE oce ! ocean dynamics and tracers 26 26 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 ! 27 31 USE in_out_manager ! I/O manager 28 32 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 133 137 INTEGER :: iif, iil, ii0, ii1, ii ! local integers 134 138 INTEGER :: ijf, ijl, ij0, ij1 ! - - 139 INTEGER :: jiseed, jjseed ! - - 135 140 INTEGER :: ios 136 141 INTEGER :: isrow ! index for ORCA1 starting row … … 187 192 END DO 188 193 189 ! (ISF) define barotropic mask and mask the ice shelf point190 ssmask(:,:)=tmask(:,:,1) ! at this stage ice shelf is not masked191 192 194 DO jk = 1, jpk 193 195 DO jj = 1, jpj … … 199 201 END DO 200 202 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 201 245 202 246 !!gm ???? -
branches/NERC/dev_r5518_GO6_r9321_isf_merge_UKESM_AIS/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90
r6486 r14573 28 28 CONTAINS 29 29 30 SUBROUTINE dom_ngb( plon, plat, kii, kjj, cdgrid )30 SUBROUTINE dom_ngb( plon, plat, kii, kjj, cdgrid, lwet, kkk ) 31 31 !!---------------------------------------------------------------------- 32 32 !! *** ROUTINE dom_ngb *** 33 33 !! 34 !! ** Purpose : find the closest grid point from a given lon/lat position34 !! ** Purpose : find the closest wet grid point from a given lon/lat position 35 35 !! 36 36 !! ** Method : look for minimum distance in cylindrical projection … … 40 40 REAL(wp) , INTENT(in ) :: plon, plat ! longitude,latitude of the point 41 41 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) 42 45 CHARACTER(len=1), INTENT(in ) :: cdgrid ! grid name 'T', 'U', 'V', 'W' 43 46 ! 47 INTEGER :: ik ! working level 44 48 INTEGER , DIMENSION(2) :: iloc 45 49 REAL(wp) :: zlon, zmini 46 50 REAL(wp), POINTER, DIMENSION(:,:) :: zglam, zgphi, zmask, zdist 51 LOGICAL :: llwet ! working logical 47 52 !!-------------------------------------------------------------------- 48 53 ! … … 51 56 CALL wrk_alloc( jpi, jpj, zglam, zgphi, zmask, zdist ) 52 57 ! 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 53 71 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 ! 67 96 zgphi(:,:) = zgphi(:,:) - plat 68 97 zdist(:,:) = zglam(:,:) * zglam(:,:) + zgphi(:,:) * zgphi(:,:) 69 98 ! 99 ! find min 70 100 IF( lk_mpp ) THEN 71 101 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 205 205 zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp ) 206 206 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 207 210 208 211 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 50 50 PUBLIC dom_zgr ! called by dom_init.F90 51 51 52 ! remove subglacial lake under the ice sheet. 53 LOGICAL , PUBLIC :: ln_isfsubgl 54 REAL(wp), PUBLIC :: rn_isfsubgllon, rn_isfsubgllat 52 55 ! !!* Namelist namzgr_sco * 53 56 LOGICAL :: ln_s_sh94 ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T) … … 395 398 IF(lwp) WRITE(numout,*) ' zgr_bat : defines level and meter bathymetry' 396 399 IF(lwp) WRITE(numout,*) ' ~~~~~~~' 400 ! 401 ! (ISF) initialisation ice shelf draft and top level 402 risfdep(:,:)=0._wp 403 misfdep(:,:)=1 397 404 ! ! ================== ! 398 405 IF( ntopo == 0 .OR. ntopo == -1 ) THEN ! defined by hand ! … … 484 491 END DO 485 492 END DO 486 risfdep(:,:)=0.e0487 misfdep(:,:)=1488 493 ! 489 494 DEALLOCATE( idta, zdta ) … … 535 540 CALL iom_close( inum ) 536 541 ! 537 risfdep(:,:)=0._wp538 misfdep(:,:)=1539 542 IF ( ln_isfcav ) THEN 540 543 CALL iom_open ( 'isf_draft_meter.nc', inum ) … … 1260 1263 !! ** Purpose : check the bathymetry in levels 1261 1264 !! 1262 !! ** Method : T He water column have to contained at least 2 cells1265 !! ** Method : The water column have to contained at least 2 cells 1263 1266 !! Bathymetry and isfdraft are modified (dig/close) to respect 1264 1267 !! 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 1266 1283 !! 1267 1284 !! ** Action : - test compatibility between isfdraft and bathy 1268 1285 !! - bathy and isfdraft are modified 1269 1286 !!---------------------------------------------------------------------- 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) 1284 1307 !!--------------------------------------------------------------------- 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 1285 1311 ! 1286 1312 IF( nn_timing == 1 ) CALL timing_start('zgr_isf') 1287 1313 ! 1288 CALL wrk_alloc( jpi, jpj, z bathy, zmask, zrisfdep)1314 CALL wrk_alloc( jpi, jpj, zdummy, zmask, zrisfdep) 1289 1315 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) 1320 901 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 ) 1324 902 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) 1298 1435 ! find the first ocean level such that the first level thickness 1299 1436 ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where 1300 1437 ! 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 ! 1301 1444 DO jk = 2, jpkm1 1302 1445 zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1303 WHERE( 0._wp < risfdep(:,:).AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+11446 WHERE( risfdep(:,:) > 0.0 .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+1 1304 1447 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 1316 1525 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) 1379 1568 END IF 1380 1569 END IF 1381 1570 END DO 1382 1571 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 1426 1575 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(:,:) ) 1430 1578 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. ) 1435 1580 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 1483 1615 ENDDO 1484 1616 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 1518 1622 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(:,:) ) 1522 1625 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 ! 1758 1633 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf') 1759 1634
Note: See TracChangeset
for help on using the changeset viewer.