- Timestamp:
- 2015-12-17T15:19:01+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_agrif_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r4982 r6092 17 17 !! 3.4 ! 2012-08 (J. Siddorn) added Siddorn and Furner stretching function 18 18 !! 3.4 ! 2012-12 (R. Bourdalle-Badie and G. Reffray) modify C1D case 19 !! 3.6 ! 2014-11 (P. Mathiot and C. Harris) add ice shelf capabilitye 19 20 !!---------------------------------------------------------------------- 20 21 … … 101 102 INTEGER :: ios 102 103 ! 103 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco 104 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav 104 105 !!---------------------------------------------------------------------- 105 106 ! … … 120 121 WRITE(numout,*) '~~~~~~~' 121 122 WRITE(numout,*) ' Namelist namzgr : set vertical coordinate' 122 WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco 123 WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps 124 WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco 123 WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco 124 WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps 125 WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco 126 WRITE(numout,*) ' ice shelf cavities ln_isfcav = ', ln_isfcav 125 127 ENDIF 126 128 … … 145 147 IF( .NOT.lk_c1d ) CALL zgr_bat_ctl ! check bathymetry (mbathy) and suppress isolated ocean points 146 148 CALL zgr_bot_level ! deepest ocean level for t-, u- and v-points 149 CALL zgr_top_level ! shallowest ocean level for T-, U-, V- points 147 150 ! 148 151 IF( lk_c1d ) THEN ! 1D config.: same mbathy value over the 3x3 domain … … 309 312 ENDIF 310 313 314 IF ( ln_isfcav ) THEN 315 ! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth) 316 ! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively 317 DO jk = 1, jpkm1 318 e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk) 319 END DO 320 e3t_1d(jpk) = e3t_1d(jpk-1) ! we don't care because this level is masked in NEMO 321 322 DO jk = 2, jpk 323 e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1) 324 END DO 325 e3w_1d(1 ) = 2._wp * (gdept_1d(1) - gdepw_1d(1)) 326 END IF 327 311 328 !!gm BUG in s-coordinate this does not work! 312 329 ! deepest/shallowest W level Above/Below ~10m … … 364 381 INTEGER :: ji, jj, jl, jk ! dummy loop indices 365 382 INTEGER :: inum ! temporary logical unit 383 INTEGER :: ierror ! error flag 366 384 INTEGER :: ii_bump, ij_bump, ih ! bump center position 367 385 INTEGER :: ii0, ii1, ij0, ij1, ik ! local indices 368 386 REAL(wp) :: r_bump , h_bump , h_oce ! bump characteristics 369 387 REAL(wp) :: zi, zj, zh, zhmin ! local scalars 370 INTEGER , POINTER, DIMENSION(:,:) :: idta ! global domain integer data371 REAL(wp), POINTER, DIMENSION(:,:) :: zdta ! global domain scalar data388 INTEGER , ALLOCATABLE, DIMENSION(:,:) :: idta ! global domain integer data 389 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdta ! global domain scalar data 372 390 !!---------------------------------------------------------------------- 373 391 ! 374 392 IF( nn_timing == 1 ) CALL timing_start('zgr_bat') 375 !376 CALL wrk_alloc( jpidta, jpjdta, idta )377 CALL wrk_alloc( jpidta, jpjdta, zdta )378 393 ! 379 394 IF(lwp) WRITE(numout,*) … … 384 399 ! ! ================== ! 385 400 ! ! global domain level and meter bathymetry (idta,zdta) 401 ! 402 ALLOCATE( idta(jpidta,jpjdta), STAT=ierror ) 403 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate idta array' ) 404 ALLOCATE( zdta(jpidta,jpjdta), STAT=ierror ) 405 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate zdta array' ) 386 406 ! 387 407 IF( ntopo == 0 ) THEN ! flat basin … … 464 484 END DO 465 485 END DO 486 risfdep(:,:)=0.e0 487 misfdep(:,:)=1 488 ! 489 DEALLOCATE( idta, zdta ) 466 490 ! 467 491 ! ! ================ ! … … 504 528 IF( ln_zps .OR. ln_sco ) THEN ! zps or sco : read meter bathymetry 505 529 CALL iom_open ( 'bathy_meter.nc', inum ) 506 CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy ) 530 IF ( ln_isfcav ) THEN 531 CALL iom_get ( inum, jpdom_data, 'Bathymetry_isf', bathy, lrowattr=.false. ) 532 ELSE 533 CALL iom_get ( inum, jpdom_data, 'Bathymetry' , bathy, lrowattr=ln_use_jattr ) 534 END IF 507 535 CALL iom_close( inum ) 508 536 ! 537 risfdep(:,:)=0._wp 538 misfdep(:,:)=1 539 IF ( ln_isfcav ) THEN 540 CALL iom_open ( 'isf_draft_meter.nc', inum ) 541 CALL iom_get ( inum, jpdom_data, 'isf_draft', risfdep ) 542 CALL iom_close( inum ) 543 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 544 END IF 545 ! 509 546 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 510 547 ! … … 544 581 ! 545 582 IF ( .not. ln_sco ) THEN !== set a minimum depth ==! 583 ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin 584 IF ( ln_isfcav ) THEN 585 WHERE (bathy == risfdep) 586 bathy = 0.0_wp ; risfdep = 0.0_wp 587 END WHERE 588 END IF 589 ! end patch 546 590 IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level 547 591 ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth … … 553 597 IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik 554 598 ENDIF 555 !556 CALL wrk_dealloc( jpidta, jpjdta, idta )557 CALL wrk_dealloc( jpidta, jpjdta, zdta )558 599 ! 559 600 IF( nn_timing == 1 ) CALL timing_stop('zgr_bat') … … 797 838 END SUBROUTINE zgr_bot_level 798 839 840 SUBROUTINE zgr_top_level 841 !!---------------------------------------------------------------------- 842 !! *** ROUTINE zgr_bot_level *** 843 !! 844 !! ** Purpose : defines the vertical index of ocean top (mik. arrays) 845 !! 846 !! ** Method : computes from misfdep with a minimum value of 1 847 !! 848 !! ** Action : mikt, miku, mikv : vertical indices of the shallowest 849 !! ocean level at t-, u- & v-points 850 !! (min value = 1) 851 !!---------------------------------------------------------------------- 852 !! 853 INTEGER :: ji, jj ! dummy loop indices 854 REAL(wp), POINTER, DIMENSION(:,:) :: zmik 855 !!---------------------------------------------------------------------- 856 ! 857 IF( nn_timing == 1 ) CALL timing_start('zgr_top_level') 858 ! 859 CALL wrk_alloc( jpi, jpj, zmik ) 860 ! 861 IF(lwp) WRITE(numout,*) 862 IF(lwp) WRITE(numout,*) ' zgr_top_level : ocean top k-index of T-, U-, V- and W-levels ' 863 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~' 864 ! 865 mikt(:,:) = MAX( misfdep(:,:) , 1 ) ! top k-index of T-level (=1) 866 ! ! top k-index of W-level (=mikt) 867 DO jj = 1, jpjm1 ! top k-index of U- (U-) level 868 DO ji = 1, jpim1 869 miku(ji,jj) = MAX( mikt(ji+1,jj ) , mikt(ji,jj) ) 870 mikv(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj) ) 871 mikf(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj), mikt(ji+1,jj ), mikt(ji+1,jj+1) ) 872 END DO 873 END DO 874 875 ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 876 zmik(:,:) = REAL( miku(:,:), wp ) ; CALL lbc_lnk(zmik,'U',1.) ; miku (:,:) = MAX( INT( zmik(:,:) ), 1 ) 877 zmik(:,:) = REAL( mikv(:,:), wp ) ; CALL lbc_lnk(zmik,'V',1.) ; mikv (:,:) = MAX( INT( zmik(:,:) ), 1 ) 878 zmik(:,:) = REAL( mikf(:,:), wp ) ; CALL lbc_lnk(zmik,'F',1.) ; mikf (:,:) = MAX( INT( zmik(:,:) ), 1 ) 879 ! 880 CALL wrk_dealloc( jpi, jpj, zmik ) 881 ! 882 IF( nn_timing == 1 ) CALL timing_stop('zgr_top_level') 883 ! 884 END SUBROUTINE zgr_top_level 799 885 800 886 SUBROUTINE zgr_zco … … 876 962 !! 877 963 INTEGER :: ji, jj, jk ! dummy loop indices 878 INTEGER :: ik, it 964 INTEGER :: ik, it, ikb, ikt ! temporary integers 879 965 LOGICAL :: ll_print ! Allow control print for debugging 880 966 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points … … 920 1006 WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1 921 1007 END DO 1008 1009 IF ( ln_isfcav ) CALL zgr_isf 922 1010 923 1011 ! Scale factors and depth at T- and W-points … … 952 1040 !gm Bug? check the gdepw_1d 953 1041 ! ... on ik 954 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0 955 & * ((gdept_1d(ik ) - gdepw_1d(ik) ) &956 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ))957 e3t_0 (ji,jj,ik) = e3t_1d(ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) &958 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )1042 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 1043 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1044 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1045 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 1046 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 959 1047 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 960 1048 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) … … 987 1075 END DO 988 1076 END DO 1077 ! 1078 IF ( ln_isfcav ) THEN 1079 ! (ISF) Definition of e3t, u, v, w for ISF case 1080 DO jj = 1, jpj 1081 DO ji = 1, jpi 1082 ik = misfdep(ji,jj) 1083 IF( ik > 1 ) THEN ! ice shelf point only 1084 IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik) 1085 gdepw_0(ji,jj,ik) = risfdep(ji,jj) 1086 !gm Bug? check the gdepw_0 1087 ! ... on ik 1088 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) & 1089 & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) & 1090 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) 1091 e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 1092 e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 1093 1094 IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column) 1095 e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 1096 ENDIF 1097 ! ... on ik / ik-1 1098 e3w_0 (ji,jj,ik ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 1099 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 1100 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 1101 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1102 ENDIF 1103 END DO 1104 END DO 1105 ! 1106 it = 0 1107 DO jj = 1, jpj 1108 DO ji = 1, jpi 1109 ik = misfdep(ji,jj) 1110 IF( ik > 1 ) THEN ! ice shelf point only 1111 e3tp (ji,jj) = e3t_0(ji,jj,ik ) 1112 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 1113 ! test 1114 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik ) 1115 IF( zdiff <= 0. .AND. lwp ) THEN 1116 it = it + 1 1117 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1118 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 1119 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1120 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj) 1121 ENDIF 1122 ENDIF 1123 END DO 1124 END DO 1125 END IF 1126 ! END (ISF) 989 1127 990 1128 ! Scale factors and depth at U-, V-, UW and VW-points … … 1005 1143 END DO 1006 1144 END DO 1145 IF ( ln_isfcav ) THEN 1146 ! (ISF) define e3uw (adapted for 2 cells in the water column) 1147 DO jj = 2, jpjm1 1148 DO ji = 2, fs_jpim1 ! vector opt. 1149 ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj)) 1150 ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj)) 1151 IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji+1,jj ,ikb ) ) & 1152 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj ,ikb-1) ) 1153 ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1)) 1154 ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1)) 1155 IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji ,jj+1,ikb ) ) & 1156 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji ,jj+1,ikb-1) ) 1157 END DO 1158 END DO 1159 END IF 1160 1007 1161 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions 1008 1162 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp ) … … 1046 1200 1047 1201 ! Compute gdep3w_0 (vertical sum of e3w) 1048 gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 1049 DO jk = 2, jpk 1050 gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk) 1051 END DO 1052 1202 IF ( ln_isfcav ) THEN ! if cavity 1203 WHERE (misfdep == 0) misfdep = 1 1204 DO jj = 1,jpj 1205 DO ji = 1,jpi 1206 gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 1207 DO jk = 2, misfdep(ji,jj) 1208 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1209 END DO 1210 IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 1211 DO jk = misfdep(ji,jj) + 1, jpk 1212 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1213 END DO 1214 END DO 1215 END DO 1216 ELSE ! no cavity 1217 gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 1218 DO jk = 2, jpk 1219 gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk) 1220 END DO 1221 END IF 1053 1222 ! ! ================= ! 1054 1223 IF(lwp .AND. ll_print) THEN ! Control print ! … … 1084 1253 ! 1085 1254 END SUBROUTINE zgr_zps 1255 1256 SUBROUTINE zgr_isf 1257 !!---------------------------------------------------------------------- 1258 !! *** ROUTINE zgr_isf *** 1259 !! 1260 !! ** Purpose : check the bathymetry in levels 1261 !! 1262 !! ** Method : THe water column have to contained at least 2 cells 1263 !! Bathymetry and isfdraft are modified (dig/close) to respect 1264 !! this criterion. 1265 !! 1266 !! 1267 !! ** Action : - test compatibility between isfdraft and bathy 1268 !! - bathy and isfdraft are modified 1269 !!---------------------------------------------------------------------- 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) 1284 !!--------------------------------------------------------------------- 1285 ! 1286 IF( nn_timing == 1 ) CALL timing_start('zgr_isf') 1287 ! 1288 CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep) 1289 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) 1298 ! find the first ocean level such that the first level thickness 1299 ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where 1300 ! e3t_0 is the reference level thickness 1301 DO jk = 2, jpkm1 1302 zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1303 WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+1 1304 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 1316 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 1379 END IF 1380 END IF 1381 END DO 1382 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 1426 IF( lk_mpp ) THEN 1427 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1428 CALL lbc_lnk( zbathy, 'T', 1. ) 1429 misfdep(:,:) = INT( zbathy(:,:) ) 1430 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(:,:) ) 1435 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 1483 ENDDO 1484 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 1518 IF( lk_mpp ) THEN 1519 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1520 CALL lbc_lnk( zbathy, 'T', 1. ) 1521 misfdep(:,:) = INT( zbathy(:,:) ) 1522 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 1758 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf') 1759 1760 END SUBROUTINE 1086 1761 1087 1762 SUBROUTINE zgr_sco
Note: See TracChangeset
for help on using the changeset viewer.
