Changeset 5837 for branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
- Timestamp:
- 2015-10-26T15:59:39+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r4624 r5837 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 … … 295 298 ENDIF 296 299 300 IF ( ln_isfcav ) THEN 301 ! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth) 302 ! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively 303 DO jk = 1, jpkm1 304 e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk) 305 END DO 306 e3t_1d(jpk) = e3t_1d(jpk-1) ! we don't care because this level is masked in NEMO 307 308 DO jk = 2, jpk 309 e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1) 310 END DO 311 e3w_1d(1 ) = 2._wp * (gdept_1d(1) - gdepw_1d(1)) 312 END IF 313 297 314 !!gm BUG in s-coordinate this does not work! 298 315 ! deepest/shallowest W level Above/Below ~10m … … 350 367 INTEGER :: ji, jj, jl, jk ! dummy loop indices 351 368 INTEGER :: inum ! temporary logical unit 369 INTEGER :: ierror ! error flag 352 370 INTEGER :: ii_bump, ij_bump, ih ! bump center position 353 371 INTEGER :: ii0, ii1, ij0, ij1, ik ! local indices 354 372 REAL(wp) :: r_bump , h_bump , h_oce ! bump characteristics 355 373 REAL(wp) :: zi, zj, zh, zhmin ! local scalars 356 INTEGER , POINTER, DIMENSION(:,:) :: idta ! global domain integer data357 REAL(wp), POINTER, DIMENSION(:,:) :: zdta ! global domain scalar data374 INTEGER , ALLOCATABLE, DIMENSION(:,:) :: idta ! global domain integer data 375 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdta ! global domain scalar data 358 376 !!---------------------------------------------------------------------- 359 377 ! 360 378 IF( nn_timing == 1 ) CALL timing_start('zgr_bat') 361 !362 CALL wrk_alloc( jpidta, jpjdta, idta )363 CALL wrk_alloc( jpidta, jpjdta, zdta )364 379 ! 365 380 IF(lwp) WRITE(numout,*) … … 370 385 ! ! ================== ! 371 386 ! ! global domain level and meter bathymetry (idta,zdta) 387 ! 388 ALLOCATE( idta(jpidta,jpjdta), STAT=ierror ) 389 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate idta array' ) 390 ALLOCATE( zdta(jpidta,jpjdta), STAT=ierror ) 391 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate zdta array' ) 372 392 ! 373 393 IF( ntopo == 0 ) THEN ! flat basin … … 450 470 END DO 451 471 END DO 472 risfdep(:,:)=0.e0 473 misfdep(:,:)=1 474 ! 475 DEALLOCATE( idta, zdta ) 452 476 ! 453 477 ! ! ================ ! … … 490 514 IF( ln_zps .OR. ln_sco ) THEN ! zps or sco : read meter bathymetry 491 515 CALL iom_open ( 'bathy_meter.nc', inum ) 492 CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy ) 516 IF ( ln_isfcav ) THEN 517 CALL iom_get ( inum, jpdom_data, 'Bathymetry_isf', bathy, lrowattr=.false. ) 518 ELSE 519 CALL iom_get ( inum, jpdom_data, 'Bathymetry' , bathy, lrowattr=ln_use_jattr ) 520 END IF 493 521 CALL iom_close( inum ) 494 522 ! 523 risfdep(:,:)=0._wp 524 misfdep(:,:)=1 525 IF ( ln_isfcav ) THEN 526 CALL iom_open ( 'isf_draft_meter.nc', inum ) 527 CALL iom_get ( inum, jpdom_data, 'isf_draft', risfdep ) 528 CALL iom_close( inum ) 529 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 530 END IF 531 ! 495 532 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 496 533 ! … … 530 567 ! 531 568 IF ( .not. ln_sco ) THEN !== set a minimum depth ==! 569 ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin 570 IF ( ln_isfcav ) THEN 571 WHERE (bathy == risfdep) 572 bathy = 0.0_wp ; risfdep = 0.0_wp 573 END WHERE 574 END IF 575 ! end patch 532 576 IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level 533 577 ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth … … 539 583 IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik 540 584 ENDIF 541 !542 CALL wrk_dealloc( jpidta, jpjdta, idta )543 CALL wrk_dealloc( jpidta, jpjdta, zdta )544 585 ! 545 586 IF( nn_timing == 1 ) CALL timing_stop('zgr_bat') … … 783 824 END SUBROUTINE zgr_bot_level 784 825 826 SUBROUTINE zgr_top_level 827 !!---------------------------------------------------------------------- 828 !! *** ROUTINE zgr_bot_level *** 829 !! 830 !! ** Purpose : defines the vertical index of ocean top (mik. arrays) 831 !! 832 !! ** Method : computes from misfdep with a minimum value of 1 833 !! 834 !! ** Action : mikt, miku, mikv : vertical indices of the shallowest 835 !! ocean level at t-, u- & v-points 836 !! (min value = 1) 837 !!---------------------------------------------------------------------- 838 !! 839 INTEGER :: ji, jj ! dummy loop indices 840 REAL(wp), POINTER, DIMENSION(:,:) :: zmik 841 !!---------------------------------------------------------------------- 842 ! 843 IF( nn_timing == 1 ) CALL timing_start('zgr_top_level') 844 ! 845 CALL wrk_alloc( jpi, jpj, zmik ) 846 ! 847 IF(lwp) WRITE(numout,*) 848 IF(lwp) WRITE(numout,*) ' zgr_top_level : ocean top k-index of T-, U-, V- and W-levels ' 849 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~' 850 ! 851 mikt(:,:) = MAX( misfdep(:,:) , 1 ) ! top k-index of T-level (=1) 852 ! ! top k-index of W-level (=mikt) 853 DO jj = 1, jpjm1 ! top k-index of U- (U-) level 854 DO ji = 1, jpim1 855 miku(ji,jj) = MAX( mikt(ji+1,jj ) , mikt(ji,jj) ) 856 mikv(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj) ) 857 mikf(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj), mikt(ji+1,jj ), mikt(ji+1,jj+1) ) 858 END DO 859 END DO 860 861 ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 862 zmik(:,:) = REAL( miku(:,:), wp ) ; CALL lbc_lnk(zmik,'U',1.) ; miku (:,:) = MAX( INT( zmik(:,:) ), 1 ) 863 zmik(:,:) = REAL( mikv(:,:), wp ) ; CALL lbc_lnk(zmik,'V',1.) ; mikv (:,:) = MAX( INT( zmik(:,:) ), 1 ) 864 zmik(:,:) = REAL( mikf(:,:), wp ) ; CALL lbc_lnk(zmik,'F',1.) ; mikf (:,:) = MAX( INT( zmik(:,:) ), 1 ) 865 ! 866 CALL wrk_dealloc( jpi, jpj, zmik ) 867 ! 868 IF( nn_timing == 1 ) CALL timing_stop('zgr_top_level') 869 ! 870 END SUBROUTINE zgr_top_level 785 871 786 872 SUBROUTINE zgr_zco … … 862 948 !! 863 949 INTEGER :: ji, jj, jk ! dummy loop indices 864 INTEGER :: ik, it 950 INTEGER :: ik, it, ikb, ikt ! temporary integers 865 951 LOGICAL :: ll_print ! Allow control print for debugging 866 952 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points … … 906 992 WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1 907 993 END DO 994 995 IF ( ln_isfcav ) CALL zgr_isf 908 996 909 997 ! Scale factors and depth at T- and W-points … … 938 1026 !gm Bug? check the gdepw_1d 939 1027 ! ... on ik 940 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0 941 & * ((gdept_1d(ik ) - gdepw_1d(ik) ) &942 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ))943 e3t_0 (ji,jj,ik) = e3t_1d(ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) &944 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )1028 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 1029 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1030 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1031 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 1032 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 945 1033 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 946 1034 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) … … 973 1061 END DO 974 1062 END DO 1063 ! 1064 IF ( ln_isfcav ) THEN 1065 ! (ISF) Definition of e3t, u, v, w for ISF case 1066 DO jj = 1, jpj 1067 DO ji = 1, jpi 1068 ik = misfdep(ji,jj) 1069 IF( ik > 1 ) THEN ! ice shelf point only 1070 IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik) 1071 gdepw_0(ji,jj,ik) = risfdep(ji,jj) 1072 !gm Bug? check the gdepw_0 1073 ! ... on ik 1074 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) & 1075 & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) & 1076 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) 1077 e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 1078 e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 1079 1080 IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column) 1081 e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 1082 ENDIF 1083 ! ... on ik / ik-1 1084 e3w_0 (ji,jj,ik ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 1085 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 1086 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 1087 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1088 ENDIF 1089 END DO 1090 END DO 1091 ! 1092 it = 0 1093 DO jj = 1, jpj 1094 DO ji = 1, jpi 1095 ik = misfdep(ji,jj) 1096 IF( ik > 1 ) THEN ! ice shelf point only 1097 e3tp (ji,jj) = e3t_0(ji,jj,ik ) 1098 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 1099 ! test 1100 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik ) 1101 IF( zdiff <= 0. .AND. lwp ) THEN 1102 it = it + 1 1103 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1104 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 1105 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1106 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj) 1107 ENDIF 1108 ENDIF 1109 END DO 1110 END DO 1111 END IF 1112 ! END (ISF) 975 1113 976 1114 ! Scale factors and depth at U-, V-, UW and VW-points … … 991 1129 END DO 992 1130 END DO 1131 IF ( ln_isfcav ) THEN 1132 ! (ISF) define e3uw (adapted for 2 cells in the water column) 1133 DO jj = 2, jpjm1 1134 DO ji = 2, fs_jpim1 ! vector opt. 1135 ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj)) 1136 ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj)) 1137 IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji+1,jj ,ikb ) ) & 1138 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj ,ikb-1) ) 1139 ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1)) 1140 ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1)) 1141 IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji ,jj+1,ikb ) ) & 1142 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji ,jj+1,ikb-1) ) 1143 END DO 1144 END DO 1145 END IF 1146 993 1147 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions 994 1148 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp ) … … 1032 1186 1033 1187 ! Compute gdep3w_0 (vertical sum of e3w) 1034 gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 1035 DO jk = 2, jpk 1036 gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk) 1037 END DO 1038 1188 IF ( ln_isfcav ) THEN ! if cavity 1189 WHERE (misfdep == 0) misfdep = 1 1190 DO jj = 1,jpj 1191 DO ji = 1,jpi 1192 gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 1193 DO jk = 2, misfdep(ji,jj) 1194 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1195 END DO 1196 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)) 1197 DO jk = misfdep(ji,jj) + 1, jpk 1198 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1199 END DO 1200 END DO 1201 END DO 1202 ELSE ! no cavity 1203 gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 1204 DO jk = 2, jpk 1205 gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk) 1206 END DO 1207 END IF 1039 1208 ! ! ================= ! 1040 1209 IF(lwp .AND. ll_print) THEN ! Control print ! … … 1070 1239 ! 1071 1240 END SUBROUTINE zgr_zps 1241 1242 SUBROUTINE zgr_isf 1243 !!---------------------------------------------------------------------- 1244 !! *** ROUTINE zgr_isf *** 1245 !! 1246 !! ** Purpose : check the bathymetry in levels 1247 !! 1248 !! ** Method : THe water column have to contained at least 2 cells 1249 !! Bathymetry and isfdraft are modified (dig/close) to respect 1250 !! this criterion. 1251 !! 1252 !! 1253 !! ** Action : - test compatibility between isfdraft and bathy 1254 !! - bathy and isfdraft are modified 1255 !!---------------------------------------------------------------------- 1256 !! 1257 INTEGER :: ji, jj, jk, jl ! dummy loop indices 1258 INTEGER :: ik, it ! temporary integers 1259 INTEGER :: id, jd, nprocd 1260 INTEGER :: icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1 ! (ISF) 1261 LOGICAL :: ll_print ! Allow control print for debugging 1262 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 1263 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t 1264 REAL(wp) :: zmax, zmin ! Maximum and minimum depth 1265 REAL(wp) :: zdiff ! temporary scalar 1266 REAL(wp) :: zrefdep ! temporary scalar 1267 REAL(wp) :: zbathydiff, zrisfdepdiff ! isf temporary scalar 1268 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH) 1269 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH) 1270 !!--------------------------------------------------------------------- 1271 ! 1272 IF( nn_timing == 1 ) CALL timing_start('zgr_isf') 1273 ! 1274 CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep) 1275 CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy ) 1276 1277 1278 ! (ISF) compute misfdep 1279 WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ; misfdep(:,:) = 1 ! open water : set misfdep to 1 1280 ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level 1281 END WHERE 1282 1283 ! Compute misfdep for ocean points (i.e. first wet level) 1284 ! find the first ocean level such that the first level thickness 1285 ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where 1286 ! e3t_0 is the reference level thickness 1287 DO jk = 2, jpkm1 1288 zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1289 WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+1 1290 END DO 1291 WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp) 1292 risfdep(:,:) = 0. ; misfdep(:,:) = 1 1293 END WHERE 1294 1295 ! 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 1296 icompt = 0 1297 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 1298 DO jl = 1, 10 1299 WHERE (bathy(:,:) .EQ. risfdep(:,:) ) 1300 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 1301 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 1302 END WHERE 1303 WHERE (mbathy(:,:) <= 0) 1304 misfdep(:,:) = 0; risfdep(:,:) = 0._wp 1305 mbathy (:,:) = 0; bathy (:,:) = 0._wp 1306 ENDWHERE 1307 IF( lk_mpp ) THEN 1308 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1309 CALL lbc_lnk( zbathy, 'T', 1. ) 1310 misfdep(:,:) = INT( zbathy(:,:) ) 1311 CALL lbc_lnk( risfdep, 'T', 1. ) 1312 CALL lbc_lnk( bathy, 'T', 1. ) 1313 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1314 CALL lbc_lnk( zbathy, 'T', 1. ) 1315 mbathy(:,:) = INT( zbathy(:,:) ) 1316 ENDIF 1317 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 1318 misfdep( 1 ,:) = misfdep(jpim1,:) ! local domain is cyclic east-west 1319 misfdep(jpi,:) = misfdep( 2 ,:) 1320 ENDIF 1321 1322 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 1323 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west 1324 mbathy(jpi,:) = mbathy( 2 ,:) 1325 ENDIF 1326 1327 ! split last cell if possible (only where water column is 2 cell or less) 1328 DO jk = jpkm1, 1, -1 1329 zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1330 WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 1331 mbathy(:,:) = jk 1332 bathy(:,:) = zmax 1333 END WHERE 1334 END DO 1335 1336 ! split top cell if possible (only where water column is 2 cell or less) 1337 DO jk = 2, jpkm1 1338 zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1339 WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy) 1340 misfdep(:,:) = jk 1341 risfdep(:,:) = zmax 1342 END WHERE 1343 END DO 1344 1345 1346 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 1347 DO jj = 1, jpj 1348 DO ji = 1, jpi 1349 ! find the minimum change option: 1350 ! test bathy 1351 IF (risfdep(ji,jj) .GT. 1) THEN 1352 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 1353 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1354 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) & 1355 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1356 1357 IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN 1358 IF (zbathydiff .LE. zrisfdepdiff) THEN 1359 bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 1360 mbathy(ji,jj)= mbathy(ji,jj) + 1 1361 ELSE 1362 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 1363 misfdep(ji,jj) = misfdep(ji,jj) - 1 1364 END IF 1365 END IF 1366 END IF 1367 END DO 1368 END DO 1369 1370 ! At least 2 levels for water thickness at T, U, and V point. 1371 DO jj = 1, jpj 1372 DO ji = 1, jpi 1373 ! find the minimum change option: 1374 ! test bathy 1375 IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 1376 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1)& 1377 & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1378 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) & 1379 & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1380 IF (zbathydiff .LE. zrisfdepdiff) THEN 1381 mbathy(ji,jj) = mbathy(ji,jj) + 1 1382 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1383 ELSE 1384 misfdep(ji,jj)= misfdep(ji,jj) - 1 1385 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 1386 END IF 1387 ENDIF 1388 END DO 1389 END DO 1390 1391 ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1) 1392 DO jj = 1, jpjm1 1393 DO ji = 1, jpim1 1394 IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 1395 zbathydiff =ABS(bathy(ji,jj ) - (gdepw_1d(mbathy (ji,jj)+1) & 1396 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat ))) 1397 zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1)) & 1398 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 1399 IF (zbathydiff .LE. zrisfdepdiff) THEN 1400 mbathy(ji,jj) = mbathy(ji,jj) + 1 1401 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) & 1402 & + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat ) 1403 ELSE 1404 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1 1405 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) & 1406 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 1407 END IF 1408 ENDIF 1409 END DO 1410 END DO 1411 1412 IF( lk_mpp ) THEN 1413 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1414 CALL lbc_lnk( zbathy, 'T', 1. ) 1415 misfdep(:,:) = INT( zbathy(:,:) ) 1416 CALL lbc_lnk( risfdep, 'T', 1. ) 1417 CALL lbc_lnk( bathy, 'T', 1. ) 1418 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1419 CALL lbc_lnk( zbathy, 'T', 1. ) 1420 mbathy(:,:) = INT( zbathy(:,:) ) 1421 ENDIF 1422 ! point V misdep(ji,jj) EQ mbathy(ji,jj+1) 1423 DO jj = 1, jpjm1 1424 DO ji = 1, jpim1 1425 IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN 1426 zbathydiff =ABS( bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) & 1427 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 1428 zrisfdepdiff=ABS(risfdep(ji,jj ) - (gdepw_1d(misfdep(ji,jj ) ) & 1429 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat ))) 1430 IF (zbathydiff .LE. zrisfdepdiff) THEN 1431 mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 1432 bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 1433 ELSE 1434 misfdep(ji,jj) = misfdep(ji,jj) - 1 1435 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat ) 1436 END IF 1437 ENDIF 1438 END DO 1439 END DO 1440 1441 1442 IF( lk_mpp ) THEN 1443 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1444 CALL lbc_lnk( zbathy, 'T', 1. ) 1445 misfdep(:,:) = INT( zbathy(:,:) ) 1446 CALL lbc_lnk( risfdep, 'T', 1. ) 1447 CALL lbc_lnk( bathy, 'T', 1. ) 1448 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1449 CALL lbc_lnk( zbathy, 'T', 1. ) 1450 mbathy(:,:) = INT( zbathy(:,:) ) 1451 ENDIF 1452 1453 ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1) 1454 DO jj = 1, jpjm1 1455 DO ji = 1, jpim1 1456 IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 1457 zbathydiff =ABS( bathy(ji ,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 1458 & + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat ))) 1459 zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) & 1460 & - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 1461 IF (zbathydiff .LE. zrisfdepdiff) THEN 1462 mbathy(ji,jj) = mbathy(ji,jj) + 1 1463 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1464 ELSE 1465 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 1466 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 1467 END IF 1468 ENDIF 1469 ENDDO 1470 ENDDO 1471 1472 IF( lk_mpp ) THEN 1473 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1474 CALL lbc_lnk( zbathy, 'T', 1. ) 1475 misfdep(:,:) = INT( zbathy(:,:) ) 1476 CALL lbc_lnk( risfdep, 'T', 1. ) 1477 CALL lbc_lnk( bathy, 'T', 1. ) 1478 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1479 CALL lbc_lnk( zbathy, 'T', 1. ) 1480 mbathy(:,:) = INT( zbathy(:,:) ) 1481 ENDIF 1482 1483 ! point U misfdep(ji,jj) EQ bathy(ji,jj+1) 1484 DO jj = 1, jpjm1 1485 DO ji = 1, jpim1 1486 IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 1487 zbathydiff =ABS( bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) & 1488 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 1489 zrisfdepdiff=ABS(risfdep(ji ,jj) - (gdepw_1d(misfdep(ji ,jj) ) & 1490 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat ))) 1491 IF (zbathydiff .LE. zrisfdepdiff) THEN 1492 mbathy(ji+1,jj) = mbathy (ji+1,jj) + 1 1493 bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) & 1494 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 1495 ELSE 1496 misfdep(ji,jj) = misfdep(ji ,jj) - 1 1497 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) & 1498 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat ) 1499 END IF 1500 ENDIF 1501 ENDDO 1502 ENDDO 1503 1504 IF( lk_mpp ) THEN 1505 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1506 CALL lbc_lnk( zbathy, 'T', 1. ) 1507 misfdep(:,:) = INT( zbathy(:,:) ) 1508 CALL lbc_lnk( risfdep, 'T', 1. ) 1509 CALL lbc_lnk( bathy, 'T', 1. ) 1510 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1511 CALL lbc_lnk( zbathy, 'T', 1. ) 1512 mbathy(:,:) = INT( zbathy(:,:) ) 1513 ENDIF 1514 END DO 1515 ! end dig bathy/ice shelf to be compatible 1516 ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness 1517 DO jl = 1,20 1518 1519 ! remove single point "bay" on isf coast line in the ice shelf draft' 1520 DO jk = 2, jpk 1521 WHERE (misfdep==0) misfdep=jpk 1522 zmask=0 1523 WHERE (misfdep .LE. jk) zmask=1 1524 DO jj = 2, jpjm1 1525 DO ji = 2, jpim1 1526 IF (misfdep(ji,jj) .EQ. jk) THEN 1527 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1528 IF (ibtest .LE. 1) THEN 1529 risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1 1530 IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) misfdep(ji,jj) = jpk 1531 END IF 1532 END IF 1533 END DO 1534 END DO 1535 END DO 1536 WHERE (misfdep==jpk) 1537 misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0. 1538 END WHERE 1539 IF( lk_mpp ) THEN 1540 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1541 CALL lbc_lnk( zbathy, 'T', 1. ) 1542 misfdep(:,:) = INT( zbathy(:,:) ) 1543 CALL lbc_lnk( risfdep, 'T', 1. ) 1544 CALL lbc_lnk( bathy, 'T', 1. ) 1545 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1546 CALL lbc_lnk( zbathy, 'T', 1. ) 1547 mbathy(:,:) = INT( zbathy(:,:) ) 1548 ENDIF 1549 1550 ! remove single point "bay" on bathy coast line beneath an ice shelf' 1551 DO jk = jpk,1,-1 1552 zmask=0 1553 WHERE (mbathy .GE. jk ) zmask=1 1554 DO jj = 2, jpjm1 1555 DO ji = 2, jpim1 1556 IF (mbathy(ji,jj) .EQ. jk .AND. misfdep(ji,jj) .GE. 2) THEN 1557 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1558 IF (ibtest .LE. 1) THEN 1559 bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1 1560 IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) mbathy(ji,jj) = 0 1561 END IF 1562 END IF 1563 END DO 1564 END DO 1565 END DO 1566 WHERE (mbathy==0) 1567 misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0. 1568 END WHERE 1569 IF( lk_mpp ) THEN 1570 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1571 CALL lbc_lnk( zbathy, 'T', 1. ) 1572 misfdep(:,:) = INT( zbathy(:,:) ) 1573 CALL lbc_lnk( risfdep, 'T', 1. ) 1574 CALL lbc_lnk( bathy, 'T', 1. ) 1575 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1576 CALL lbc_lnk( zbathy, 'T', 1. ) 1577 mbathy(:,:) = INT( zbathy(:,:) ) 1578 ENDIF 1579 1580 ! fill hole in ice shelf 1581 zmisfdep = misfdep 1582 zrisfdep = risfdep 1583 WHERE (zmisfdep .LE. 1) zmisfdep=jpk 1584 DO jj = 2, jpjm1 1585 DO ji = 2, jpim1 1586 ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj ) 1587 ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1) 1588 IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj ) - 1) 1589 IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj ) - 1) 1590 IF( zmisfdep(ji,jj) .GE. mbathy(ji ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji ,jj-1) - 1) 1591 IF( zmisfdep(ji,jj) .GE. mbathy(ji ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji ,jj+1) - 1) 1592 ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1593 IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN 1594 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp 1595 END IF 1596 IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) .GE. 2) THEN 1597 misfdep(ji,jj) = ibtest 1598 risfdep(ji,jj) = gdepw_1d(ibtest) 1599 ENDIF 1600 ENDDO 1601 ENDDO 1602 1603 IF( lk_mpp ) THEN 1604 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1605 CALL lbc_lnk( zbathy, 'T', 1. ) 1606 misfdep(:,:) = INT( zbathy(:,:) ) 1607 CALL lbc_lnk( risfdep, 'T', 1. ) 1608 CALL lbc_lnk( bathy, 'T', 1. ) 1609 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1610 CALL lbc_lnk( zbathy, 'T', 1. ) 1611 mbathy(:,:) = INT( zbathy(:,:) ) 1612 ENDIF 1613 ! 1614 !! fill hole in bathymetry 1615 zmbathy (:,:)=mbathy (:,:) 1616 DO jj = 2, jpjm1 1617 DO ji = 2, jpim1 1618 ibtestim1 = zmbathy(ji-1,jj ) ; ibtestip1 = zmbathy(ji+1,jj ) 1619 ibtestjm1 = zmbathy(ji ,jj-1) ; ibtestjp1 = zmbathy(ji ,jj+1) 1620 IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj ) + 1) 1621 IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj ) ) ibtestip1 = 0 1622 IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj-1) ) ibtestjm1 = 0 1623 IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj+1) ) ibtestjp1 = 0 1624 ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1625 IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN 1626 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 1627 END IF 1628 IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) .GE. 2) THEN 1629 mbathy(ji,jj) = ibtest 1630 bathy(ji,jj) = gdepw_1d(ibtest+1) 1631 ENDIF 1632 END DO 1633 END DO 1634 IF( lk_mpp ) THEN 1635 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1636 CALL lbc_lnk( zbathy, 'T', 1. ) 1637 misfdep(:,:) = INT( zbathy(:,:) ) 1638 CALL lbc_lnk( risfdep, 'T', 1. ) 1639 CALL lbc_lnk( bathy, 'T', 1. ) 1640 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1641 CALL lbc_lnk( zbathy, 'T', 1. ) 1642 mbathy(:,:) = INT( zbathy(:,:) ) 1643 ENDIF 1644 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 1645 DO jj = 1, jpjm1 1646 DO ji = 1, jpim1 1647 IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN 1648 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1649 END IF 1650 END DO 1651 END DO 1652 IF( lk_mpp ) THEN 1653 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1654 CALL lbc_lnk( zbathy, 'T', 1. ) 1655 misfdep(:,:) = INT( zbathy(:,:) ) 1656 CALL lbc_lnk( risfdep, 'T', 1. ) 1657 CALL lbc_lnk( bathy, 'T', 1. ) 1658 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1659 CALL lbc_lnk( zbathy, 'T', 1. ) 1660 mbathy(:,:) = INT( zbathy(:,:) ) 1661 ENDIF 1662 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 1663 DO jj = 1, jpjm1 1664 DO ji = 1, jpim1 1665 IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN 1666 mbathy(ji+1,jj) = mbathy(ji+1,jj) - 1; bathy(ji+1,jj) = gdepw_1d(mbathy(ji+1,jj)+1) ; 1667 END IF 1668 END DO 1669 END DO 1670 IF( lk_mpp ) THEN 1671 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1672 CALL lbc_lnk( zbathy, 'T', 1. ) 1673 misfdep(:,:) = INT( zbathy(:,:) ) 1674 CALL lbc_lnk( risfdep, 'T', 1. ) 1675 CALL lbc_lnk( bathy, 'T', 1. ) 1676 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1677 CALL lbc_lnk( zbathy, 'T', 1. ) 1678 mbathy(:,:) = INT( zbathy(:,:) ) 1679 ENDIF 1680 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 1681 DO jj = 1, jpjm1 1682 DO ji = 1, jpi 1683 IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN 1684 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1685 END IF 1686 END DO 1687 END DO 1688 IF( lk_mpp ) THEN 1689 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1690 CALL lbc_lnk( zbathy, 'T', 1. ) 1691 misfdep(:,:) = INT( zbathy(:,:) ) 1692 CALL lbc_lnk( risfdep, 'T', 1. ) 1693 CALL lbc_lnk( bathy, 'T', 1. ) 1694 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1695 CALL lbc_lnk( zbathy, 'T', 1. ) 1696 mbathy(:,:) = INT( zbathy(:,:) ) 1697 ENDIF 1698 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 1699 DO jj = 1, jpjm1 1700 DO ji = 1, jpi 1701 IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN 1702 mbathy(ji,jj+1) = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ; 1703 END IF 1704 END DO 1705 END DO 1706 IF( lk_mpp ) THEN 1707 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1708 CALL lbc_lnk( zbathy, 'T', 1. ) 1709 misfdep(:,:) = INT( zbathy(:,:) ) 1710 CALL lbc_lnk( risfdep, 'T', 1. ) 1711 CALL lbc_lnk( bathy, 'T', 1. ) 1712 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1713 CALL lbc_lnk( zbathy, 'T', 1. ) 1714 mbathy(:,:) = INT( zbathy(:,:) ) 1715 ENDIF 1716 ! if not compatible after all check, mask T 1717 DO jj = 1, jpj 1718 DO ji = 1, jpi 1719 IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN 1720 misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0._wp ; 1721 END IF 1722 END DO 1723 END DO 1724 1725 WHERE (mbathy(:,:) == 1) 1726 mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp 1727 END WHERE 1728 END DO 1729 ! end check compatibility ice shelf/bathy 1730 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1731 WHERE (misfdep(:,:) <= 5) 1732 misfdep = 1; risfdep = 0.0_wp; 1733 END WHERE 1734 1735 IF( icompt == 0 ) THEN 1736 IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry' 1737 ELSE 1738 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry' 1739 ENDIF 1740 1741 CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 1742 CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 1743 1744 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf') 1745 1746 END SUBROUTINE 1072 1747 1073 1748 SUBROUTINE zgr_sco … … 1445 2120 DO jk = 1, jpkm1 1446 2121 IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 1447 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 01448 END DO2122 END DO 2123 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 1449 2124 END DO 1450 2125 END DO
Note: See TracChangeset
for help on using the changeset viewer.