- Timestamp:
- 2015-12-16T16:44:35+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r6060 r6069 382 382 !! - bathy : meter bathymetry (in meters) 383 383 !!---------------------------------------------------------------------- 384 INTEGER :: ji, jj, j l, jk ! dummy loop indices384 INTEGER :: ji, jj, jk ! dummy loop indices 385 385 INTEGER :: inum ! temporary logical unit 386 386 INTEGER :: ierror ! error flag … … 544 544 CALL iom_close( inum ) 545 545 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 546 547 ! set grounded point to 0 548 ! (a treshold could be set here if needed, or set it offline based on the grounded fraction) 549 WHERE ( bathy(:,:) <= risfdep(:,:) + rn_isfhmin ) 550 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 551 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 552 END WHERE 546 553 END IF 547 554 ! … … 581 588 ! 582 589 IF ( .not. ln_sco ) THEN !== set a minimum depth ==! 583 ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin584 IF ( ln_isfcav ) THEN585 WHERE (bathy == risfdep)586 bathy = 0.0_wp ; risfdep = 0.0_wp587 END WHERE588 END IF589 ! end patch590 590 IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level 591 591 ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth … … 830 830 SUBROUTINE zgr_top_level 831 831 !!---------------------------------------------------------------------- 832 !! *** ROUTINE zgr_ bot_level ***832 !! *** ROUTINE zgr_top_level *** 833 833 !! 834 834 !! ** Purpose : defines the vertical index of ocean top (mik. arrays) … … 954 954 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 955 955 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t 956 REAL(wp) :: zmax ! Maximum depth957 956 REAL(wp) :: zdiff ! temporary scalar 958 REAL(wp) :: z refdep! temporary scalar957 REAL(wp) :: zmax ! temporary scalar 959 958 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt 960 959 !!--------------------------------------------------------------------- … … 986 985 END DO 987 986 988 IF ( ln_isfcav ) CALL zgr_isf989 990 987 ! Scale factors and depth at T- and W-points 991 988 DO jk = 1, jpk ! intitialization to the reference z-coordinate … … 995 992 e3w_0 (:,:,jk) = e3w_1d (jk) 996 993 END DO 994 995 ! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf 996 IF ( ln_isfcav ) CALL zgr_isf 997 998 ! Scale factors and depth at T- and W-points 999 IF ( .NOT. ln_isfcav ) THEN 1000 DO jj = 1, jpj 1001 DO ji = 1, jpi 1002 ik = mbathy(ji,jj) 1003 IF( ik > 0 ) THEN ! ocean point only 1004 ! max ocean level case 1005 IF( ik == jpkm1 ) THEN 1006 zdepwp = bathy(ji,jj) 1007 ze3tp = bathy(ji,jj) - gdepw_1d(ik) 1008 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 1009 e3t_0(ji,jj,ik ) = ze3tp 1010 e3t_0(ji,jj,ik+1) = ze3tp 1011 e3w_0(ji,jj,ik ) = ze3wp 1012 e3w_0(ji,jj,ik+1) = ze3tp 1013 gdepw_0(ji,jj,ik+1) = zdepwp 1014 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 1015 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 1016 ! 1017 ELSE ! standard case 1018 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 1019 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1020 ENDIF 1021 !gm Bug? check the gdepw_1d 1022 ! ... on ik 1023 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 1024 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1025 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1026 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 1027 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 1028 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 1029 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 1030 ! ... on ik+1 1031 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1032 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1033 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 1034 ENDIF 1035 ENDIF 1036 END DO 1037 END DO 1038 ! 1039 it = 0 1040 DO jj = 1, jpj 1041 DO ji = 1, jpi 1042 ik = mbathy(ji,jj) 1043 IF( ik > 0 ) THEN ! ocean point only 1044 e3tp (ji,jj) = e3t_0(ji,jj,ik) 1045 e3wp (ji,jj) = e3w_0(ji,jj,ik) 1046 ! test 1047 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik ) 1048 IF( zdiff <= 0._wp .AND. lwp ) THEN 1049 it = it + 1 1050 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1051 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 1052 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1053 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik ) 1054 ENDIF 1055 ENDIF 1056 END DO 1057 END DO 1058 END IF 1059 ! 1060 ! Scale factors and depth at U-, V-, UW and VW-points 1061 DO jk = 1, jpk ! initialisation to z-scale factors 1062 e3u_0 (:,:,jk) = e3t_1d(jk) 1063 e3v_0 (:,:,jk) = e3t_1d(jk) 1064 e3uw_0(:,:,jk) = e3w_1d(jk) 1065 e3vw_0(:,:,jk) = e3w_1d(jk) 1066 END DO 1067 1068 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors 1069 DO jj = 1, jpjm1 1070 DO ji = 1, fs_jpim1 ! vector opt. 1071 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 1072 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 1073 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 1074 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 1075 END DO 1076 END DO 1077 END DO 1078 IF ( ln_isfcav ) THEN 1079 ! (ISF) define e3uw (adapted for 2 cells in the water column) 1080 DO jj = 2, jpjm1 1081 DO ji = 2, fs_jpim1 ! vector opt. 1082 ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj)) 1083 ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj)) 1084 IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji+1,jj ,ikb ) ) & 1085 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj ,ikb-1) ) 1086 ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1)) 1087 ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1)) 1088 IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji ,jj+1,ikb ) ) & 1089 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji ,jj+1,ikb-1) ) 1090 END DO 1091 END DO 1092 END IF 1093 1094 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions 1095 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 1096 ! 1097 1098 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1099 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk) 1100 WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk) 1101 WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk) 1102 WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk) 1103 END DO 1104 1105 ! Scale factor at F-point 1106 DO jk = 1, jpk ! initialisation to z-scale factors 1107 e3f_0(:,:,jk) = e3t_1d(jk) 1108 END DO 1109 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors 1110 DO jj = 1, jpjm1 1111 DO ji = 1, fs_jpim1 ! vector opt. 1112 e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 1113 END DO 1114 END DO 1115 END DO 1116 CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions 1117 ! 1118 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1119 WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk) 1120 END DO 1121 !!gm bug ? : must be a do loop with mj0,mj1 997 1122 ! 1123 e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2 1124 e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 1125 e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 1126 e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 1127 e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 1128 1129 ! Control of the sign 1130 IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' ) 1131 IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' ) 1132 IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' ) 1133 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' ) 1134 1135 ! Compute gde3w_0 (vertical sum of e3w) 1136 IF ( ln_isfcav ) THEN ! if cavity 1137 WHERE( misfdep == 0 ) misfdep = 1 1138 DO jj = 1,jpj 1139 DO ji = 1,jpi 1140 gde3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 1141 DO jk = 2, misfdep(ji,jj) 1142 gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1143 END DO 1144 IF( misfdep(ji,jj) >= 2 ) gde3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 1145 DO jk = misfdep(ji,jj) + 1, jpk 1146 gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1147 END DO 1148 END DO 1149 END DO 1150 ELSE ! no cavity 1151 gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 1152 DO jk = 2, jpk 1153 gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 1154 END DO 1155 END IF 1156 ! 1157 CALL wrk_dealloc( jpi,jpj,jpk, zprt ) 1158 ! 1159 IF( nn_timing == 1 ) CALL timing_stop('zgr_zps') 1160 ! 1161 END SUBROUTINE zgr_zps 1162 1163 1164 SUBROUTINE zgr_isf 1165 !!---------------------------------------------------------------------- 1166 !! *** ROUTINE zgr_isf *** 1167 !! 1168 !! ** Purpose : check the bathymetry in levels 1169 !! 1170 !! ** Method : THe water column have to contained at least 2 cells 1171 !! Bathymetry and isfdraft are modified (dig/close) to respect 1172 !! this criterion. 1173 !! 1174 !! ** Action : - test compatibility between isfdraft and bathy 1175 !! - bathy and isfdraft are modified 1176 !!---------------------------------------------------------------------- 1177 INTEGER :: ji, jj, jl, jk ! dummy loop indices 1178 INTEGER :: ik, it ! temporary integers 1179 INTEGER :: icompt, ibtest ! (ISF) 1180 INTEGER :: ibtestim1, ibtestip1 ! (ISF) 1181 INTEGER :: ibtestjm1, ibtestjp1 ! (ISF) 1182 REAL(wp) :: zdepth ! Ajusted ocean depth to avoid too small e3t 1183 REAL(wp) :: zmax ! Maximum and minimum depth 1184 REAL(wp) :: zbathydiff ! isf temporary scalar 1185 REAL(wp) :: zrisfdepdiff ! isf temporary scalar 1186 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 1187 REAL(wp) :: zdepwp ! Ajusted ocean depth to avoid too small e3t 1188 REAL(wp) :: zdiff ! temporary scalar 1189 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH) 1190 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH) 1191 !!--------------------------------------------------------------------- 1192 ! 1193 IF( nn_timing == 1 ) CALL timing_start('zgr_isf') 1194 ! 1195 CALL wrk_alloc( jpi,jpj, zbathy, zmask, zrisfdep) 1196 CALL wrk_alloc( jpi,jpj, zmisfdep, zmbathy ) 1197 1198 ! (ISF) compute misfdep 1199 WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) /= 0 ) ; misfdep(:,:) = 1 ! open water : set misfdep to 1 1200 ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level 1201 END WHERE 1202 1203 ! Compute misfdep for ocean points (i.e. first wet level) 1204 ! find the first ocean level such that the first level thickness 1205 ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where 1206 ! e3t_0 is the reference level thickness 1207 DO jk = 2, jpkm1 1208 zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1209 WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+1 1210 END DO 1211 WHERE ( 0._wp < risfdep(:,:) .AND. risfdep(:,:) <= e3t_1d(1) ) 1212 risfdep(:,:) = 0. ; misfdep(:,:) = 1 1213 END WHERE 1214 1215 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1216 WHERE (risfdep(:,:) <= 10._wp .AND. misfdep(:,:) > 1) 1217 misfdep = 0; risfdep = 0.0_wp; 1218 mbathy = 0; bathy = 0.0_wp; 1219 END WHERE 1220 WHERE (bathy(:,:) <= 30.0_wp .AND. gphit < -60._wp) 1221 misfdep = 0; risfdep = 0.0_wp; 1222 mbathy = 0; bathy = 0.0_wp; 1223 END WHERE 1224 1225 ! 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 1226 icompt = 0 1227 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 1228 DO jl = 1, 10 1229 ! check at each iteration if isf is grounded or not (1cm treshold have to be update after first coupling experiments) 1230 WHERE (bathy(:,:) <= risfdep(:,:) + rn_isfhmin) 1231 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 1232 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 1233 END WHERE 1234 WHERE (mbathy(:,:) <= 0) 1235 misfdep(:,:) = 0; risfdep(:,:) = 0._wp 1236 mbathy (:,:) = 0; bathy (:,:) = 0._wp 1237 END WHERE 1238 IF( lk_mpp ) THEN 1239 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1240 CALL lbc_lnk( zbathy, 'T', 1. ) 1241 misfdep(:,:) = INT( zbathy(:,:) ) 1242 1243 CALL lbc_lnk( risfdep,'T', 1. ) 1244 CALL lbc_lnk( bathy, 'T', 1. ) 1245 1246 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1247 CALL lbc_lnk( zbathy, 'T', 1. ) 1248 mbathy(:,:) = INT( zbathy(:,:) ) 1249 ENDIF 1250 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 1251 misfdep( 1 ,:) = misfdep(jpim1,:) ! local domain is cyclic east-west 1252 misfdep(jpi,:) = misfdep( 2 ,:) 1253 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west 1254 mbathy(jpi,:) = mbathy( 2 ,:) 1255 ENDIF 1256 1257 ! split last cell if possible (only where water column is 2 cell or less) 1258 ! if coupled to ice sheet, we do not modify the bathymetry (can be discuss). 1259 IF ( .NOT. ln_iscpl) THEN 1260 DO jk = jpkm1, 1, -1 1261 zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1262 WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 1263 mbathy(:,:) = jk 1264 bathy(:,:) = zmax 1265 END WHERE 1266 END DO 1267 END IF 1268 1269 ! split top cell if possible (only where water column is 2 cell or less) 1270 DO jk = 2, jpkm1 1271 zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1272 WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy) 1273 misfdep(:,:) = jk 1274 risfdep(:,:) = zmax 1275 END WHERE 1276 END DO 1277 1278 1279 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 1280 DO jj = 1, jpj 1281 DO ji = 1, jpi 1282 ! find the minimum change option: 1283 ! test bathy 1284 IF (risfdep(ji,jj) > 1) THEN 1285 IF ( .NOT. ln_iscpl ) THEN 1286 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 1287 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1288 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) & 1289 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1290 IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN 1291 IF (zbathydiff <= zrisfdepdiff) THEN 1292 bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 1293 mbathy(ji,jj)= mbathy(ji,jj) + 1 1294 ELSE 1295 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 1296 misfdep(ji,jj) = misfdep(ji,jj) - 1 1297 END IF 1298 ENDIF 1299 ELSE 1300 IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN 1301 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 1302 misfdep(ji,jj) = misfdep(ji,jj) - 1 1303 END IF 1304 END IF 1305 END IF 1306 END DO 1307 END DO 1308 1309 ! At least 2 levels for water thickness at T, U, and V point. 1310 DO jj = 1, jpj 1311 DO ji = 1, jpi 1312 ! find the minimum change option: 1313 ! test bathy 1314 IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1315 IF ( .NOT. ln_iscpl ) THEN 1316 zbathydiff =ABS(bathy(ji,jj) - ( gdepw_1d(mbathy (ji,jj)+1) & 1317 & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1318 zrisfdepdiff=ABS(risfdep(ji,jj) - ( gdepw_1d(misfdep(ji,jj) ) & 1319 & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1320 IF (zbathydiff <= zrisfdepdiff) THEN 1321 mbathy(ji,jj) = mbathy(ji,jj) + 1 1322 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1323 ELSE 1324 misfdep(ji,jj)= misfdep(ji,jj) - 1 1325 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 1326 END IF 1327 ELSE 1328 misfdep(ji,jj)= misfdep(ji,jj) - 1 1329 risfdep(ji,jj)= gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 1330 END IF 1331 ENDIF 1332 END DO 1333 END DO 1334 1335 ! point V mbathy(ji,jj) == misfdep(ji,jj+1) 1336 DO jj = 1, jpjm1 1337 DO ji = 1, jpim1 1338 IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1339 IF ( .NOT. ln_iscpl ) THEN 1340 zbathydiff =ABS(bathy(ji,jj ) - ( gdepw_1d(mbathy (ji,jj)+1) & 1341 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat ))) 1342 zrisfdepdiff=ABS(risfdep(ji,jj+1) - ( gdepw_1d(misfdep(ji,jj+1)) & 1343 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 1344 IF (zbathydiff <= zrisfdepdiff) THEN 1345 mbathy(ji,jj) = mbathy(ji,jj) + 1 1346 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat ) 1347 ELSE 1348 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1 1349 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 1350 END IF 1351 ELSE 1352 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1 1353 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 1354 END IF 1355 ENDIF 1356 END DO 1357 END DO 1358 1359 IF( lk_mpp ) THEN 1360 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1361 CALL lbc_lnk( zbathy, 'T', 1. ) 1362 misfdep(:,:) = INT( zbathy(:,:) ) 1363 1364 CALL lbc_lnk( risfdep,'T', 1. ) 1365 CALL lbc_lnk( bathy, 'T', 1. ) 1366 1367 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1368 CALL lbc_lnk( zbathy, 'T', 1. ) 1369 mbathy(:,:) = INT( zbathy(:,:) ) 1370 ENDIF 1371 ! point V misdep(ji,jj) == mbathy(ji,jj+1) 1372 DO jj = 1, jpjm1 1373 DO ji = 1, jpim1 1374 IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) > 1) THEN 1375 IF ( .NOT. ln_iscpl ) THEN 1376 zbathydiff =ABS( bathy(ji,jj+1) - ( gdepw_1d(mbathy (ji,jj+1)+1) & 1377 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+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 <= zrisfdepdiff) THEN 1381 mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 1382 bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+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 ELSE 1388 misfdep(ji,jj) = misfdep(ji,jj) - 1 1389 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat ) 1390 END IF 1391 ENDIF 1392 END DO 1393 END DO 1394 1395 1396 IF( lk_mpp ) THEN 1397 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1398 CALL lbc_lnk( zbathy, 'T', 1. ) 1399 misfdep(:,:) = INT( zbathy(:,:) ) 1400 1401 CALL lbc_lnk( risfdep,'T', 1. ) 1402 CALL lbc_lnk( bathy, 'T', 1. ) 1403 1404 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1405 CALL lbc_lnk( zbathy, 'T', 1. ) 1406 mbathy(:,:) = INT( zbathy(:,:) ) 1407 ENDIF 1408 1409 ! point U mbathy(ji,jj) == misfdep(ji,jj+1) 1410 DO jj = 1, jpjm1 1411 DO ji = 1, jpim1 1412 IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1413 IF ( .NOT. ln_iscpl ) THEN 1414 zbathydiff =ABS( bathy(ji ,jj) - ( gdepw_1d(mbathy (ji,jj)+1) & 1415 & + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat ))) 1416 zrisfdepdiff=ABS(risfdep(ji+1,jj) - ( gdepw_1d(misfdep(ji+1,jj)) & 1417 & - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 1418 IF (zbathydiff <= zrisfdepdiff) THEN 1419 mbathy(ji,jj) = mbathy(ji,jj) + 1 1420 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1421 ELSE 1422 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 1423 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 1424 END IF 1425 ELSE 1426 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 1427 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 1428 ENDIF 1429 ENDIF 1430 ENDDO 1431 ENDDO 1432 1433 IF( lk_mpp ) THEN 1434 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1435 CALL lbc_lnk( zbathy, 'T', 1. ) 1436 misfdep(:,:) = INT( zbathy(:,:) ) 1437 1438 CALL lbc_lnk( risfdep,'T', 1. ) 1439 CALL lbc_lnk( bathy, 'T', 1. ) 1440 1441 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1442 CALL lbc_lnk( zbathy, 'T', 1. ) 1443 mbathy(:,:) = INT( zbathy(:,:) ) 1444 ENDIF 1445 1446 ! point U misfdep(ji,jj) == bathy(ji,jj+1) 1447 DO jj = 1, jpjm1 1448 DO ji = 1, jpim1 1449 IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) > 1) THEN 1450 IF ( .NOT. ln_iscpl ) THEN 1451 zbathydiff =ABS( bathy(ji+1,jj) - ( gdepw_1d(mbathy (ji+1,jj)+1) & 1452 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 1453 zrisfdepdiff=ABS(risfdep(ji ,jj) - ( gdepw_1d(misfdep(ji ,jj) ) & 1454 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat ))) 1455 IF (zbathydiff <= zrisfdepdiff) THEN 1456 mbathy(ji+1,jj) = mbathy (ji+1,jj) + 1 1457 bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 1458 ELSE 1459 misfdep(ji,jj) = misfdep(ji ,jj) - 1 1460 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat ) 1461 END IF 1462 ELSE 1463 misfdep(ji,jj) = misfdep(ji ,jj) - 1 1464 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat ) 1465 ENDIF 1466 ENDIF 1467 ENDDO 1468 ENDDO 1469 1470 IF( lk_mpp ) THEN 1471 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1472 CALL lbc_lnk( zbathy, 'T', 1. ) 1473 misfdep(:,:) = INT( zbathy(:,:) ) 1474 1475 CALL lbc_lnk( risfdep,'T', 1. ) 1476 CALL lbc_lnk( bathy, 'T', 1. ) 1477 1478 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1479 CALL lbc_lnk( zbathy, 'T', 1. ) 1480 mbathy(:,:) = INT( zbathy(:,:) ) 1481 ENDIF 1482 END DO 1483 ! end dig bathy/ice shelf to be compatible 1484 ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness 1485 DO jl = 1,20 1486 1487 ! remove single point "bay" on isf coast line in the ice shelf draft' 1488 DO jk = 2, jpk 1489 WHERE (misfdep==0) misfdep=jpk 1490 zmask=0._wp 1491 WHERE (misfdep <= jk) zmask=1 1492 DO jj = 2, jpjm1 1493 DO ji = 2, jpim1 1494 IF (misfdep(ji,jj) == jk) THEN 1495 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1496 IF (ibtest <= 1) THEN 1497 risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1 1498 IF (misfdep(ji,jj) > mbathy(ji,jj)) misfdep(ji,jj) = jpk 1499 END IF 1500 END IF 1501 END DO 1502 END DO 1503 END DO 1504 WHERE (misfdep==jpk) 1505 misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 1506 END WHERE 1507 IF( lk_mpp ) THEN 1508 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1509 CALL lbc_lnk( zbathy, 'T', 1. ) 1510 misfdep(:,:) = INT( zbathy(:,:) ) 1511 1512 CALL lbc_lnk( risfdep,'T', 1. ) 1513 CALL lbc_lnk( bathy, 'T', 1. ) 1514 1515 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1516 CALL lbc_lnk( zbathy, 'T', 1. ) 1517 mbathy(:,:) = INT( zbathy(:,:) ) 1518 ENDIF 1519 1520 ! remove single point "bay" on bathy coast line beneath an ice shelf' 1521 DO jk = jpk,1,-1 1522 zmask=0._wp 1523 WHERE (mbathy >= jk ) zmask=1 1524 DO jj = 2, jpjm1 1525 DO ji = 2, jpim1 1526 IF (mbathy(ji,jj) == jk .AND. misfdep(ji,jj) >= 2) THEN 1527 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1528 IF (ibtest <= 1) THEN 1529 bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1 1530 IF (misfdep(ji,jj) > mbathy(ji,jj)) mbathy(ji,jj) = 0 1531 END IF 1532 END IF 1533 END DO 1534 END DO 1535 END DO 1536 WHERE (mbathy==0) 1537 misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 1538 END WHERE 1539 IF( lk_mpp ) THEN 1540 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1541 CALL lbc_lnk( zbathy, 'T', 1. ) 1542 misfdep(:,:) = INT( zbathy(:,:) ) 1543 1544 CALL lbc_lnk( risfdep,'T', 1. ) 1545 CALL lbc_lnk( bathy, 'T', 1. ) 1546 1547 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1548 CALL lbc_lnk( zbathy, 'T', 1. ) 1549 mbathy(:,:) = INT( zbathy(:,:) ) 1550 ENDIF 1551 1552 ! fill hole in ice shelf 1553 zmisfdep = misfdep 1554 zrisfdep = risfdep 1555 WHERE (zmisfdep <= 1._wp) zmisfdep=jpk 1556 DO jj = 2, jpjm1 1557 DO ji = 2, jpim1 1558 ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj ) 1559 ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1) 1560 IF( zmisfdep(ji,jj) >= mbathy(ji-1,jj ) ) ibtestim1 = jpk 1561 IF( zmisfdep(ji,jj) >= mbathy(ji+1,jj ) ) ibtestip1 = jpk 1562 IF( zmisfdep(ji,jj) >= mbathy(ji ,jj-1) ) ibtestjm1 = jpk 1563 IF( zmisfdep(ji,jj) >= mbathy(ji ,jj+1) ) ibtestjp1 = jpk 1564 ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1565 IF( ibtest == jpk .AND. misfdep(ji,jj) >= 2) THEN 1566 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp 1567 END IF 1568 IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) >= 2) THEN 1569 misfdep(ji,jj) = ibtest 1570 risfdep(ji,jj) = gdepw_1d(ibtest) 1571 ENDIF 1572 ENDDO 1573 ENDDO 1574 1575 IF( lk_mpp ) THEN 1576 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1577 CALL lbc_lnk( zbathy, 'T', 1. ) 1578 misfdep(:,:) = INT( zbathy(:,:) ) 1579 1580 CALL lbc_lnk( risfdep, 'T', 1. ) 1581 CALL lbc_lnk( bathy, 'T', 1. ) 1582 1583 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1584 CALL lbc_lnk( zbathy, 'T', 1. ) 1585 mbathy(:,:) = INT( zbathy(:,:) ) 1586 ENDIF 1587 ! 1588 !! fill hole in bathymetry 1589 zmbathy (:,:)=mbathy (:,:) 1590 DO jj = 2, jpjm1 1591 DO ji = 2, jpim1 1592 ibtestim1 = zmbathy(ji-1,jj ) ; ibtestip1 = zmbathy(ji+1,jj ) 1593 ibtestjm1 = zmbathy(ji ,jj-1) ; ibtestjp1 = zmbathy(ji ,jj+1) 1594 IF( zmbathy(ji,jj) < misfdep(ji-1,jj ) ) ibtestim1 = 0 1595 IF( zmbathy(ji,jj) < misfdep(ji+1,jj ) ) ibtestip1 = 0 1596 IF( zmbathy(ji,jj) < misfdep(ji ,jj-1) ) ibtestjm1 = 0 1597 IF( zmbathy(ji,jj) < misfdep(ji ,jj+1) ) ibtestjp1 = 0 1598 ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1599 IF( ibtest == 0 .AND. misfdep(ji,jj) >= 2) THEN 1600 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 1601 END IF 1602 IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) >= 2) THEN 1603 mbathy(ji,jj) = ibtest 1604 bathy(ji,jj) = gdepw_1d(ibtest+1) 1605 ENDIF 1606 END DO 1607 END DO 1608 IF( lk_mpp ) THEN 1609 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1610 CALL lbc_lnk( zbathy, 'T', 1. ) 1611 misfdep(:,:) = INT( zbathy(:,:) ) 1612 1613 CALL lbc_lnk( risfdep, 'T', 1. ) 1614 CALL lbc_lnk( bathy, 'T', 1. ) 1615 1616 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1617 CALL lbc_lnk( zbathy, 'T', 1. ) 1618 mbathy(:,:) = INT( zbathy(:,:) ) 1619 ENDIF 1620 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 1621 DO jj = 1, jpjm1 1622 DO ji = 1, jpim1 1623 IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 1624 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1625 END IF 1626 END DO 1627 END DO 1628 IF( lk_mpp ) THEN 1629 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1630 CALL lbc_lnk( zbathy, 'T', 1. ) 1631 misfdep(:,:) = INT( zbathy(:,:) ) 1632 1633 CALL lbc_lnk( risfdep, 'T', 1. ) 1634 CALL lbc_lnk( bathy, 'T', 1. ) 1635 1636 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1637 CALL lbc_lnk( zbathy, 'T', 1. ) 1638 mbathy(:,:) = INT( zbathy(:,:) ) 1639 ENDIF 1640 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 1641 DO jj = 1, jpjm1 1642 DO ji = 1, jpim1 1643 IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 1644 mbathy(ji+1,jj) = mbathy(ji+1,jj) - 1; bathy(ji+1,jj) = gdepw_1d(mbathy(ji+1,jj)+1) ; 1645 END IF 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 1653 CALL lbc_lnk( risfdep,'T', 1. ) 1654 CALL lbc_lnk( bathy, 'T', 1. ) 1655 1656 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1657 CALL lbc_lnk( zbathy, 'T', 1. ) 1658 mbathy(:,:) = INT( zbathy(:,:) ) 1659 ENDIF 1660 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 1661 DO jj = 1, jpjm1 1662 DO ji = 1, jpi 1663 IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 1664 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1665 END IF 1666 END DO 1667 END DO 1668 IF( lk_mpp ) THEN 1669 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1670 CALL lbc_lnk( zbathy, 'T', 1. ) 1671 misfdep(:,:) = INT( zbathy(:,:) ) 1672 1673 CALL lbc_lnk( risfdep,'T', 1. ) 1674 CALL lbc_lnk( bathy, 'T', 1. ) 1675 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 (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 1684 mbathy(ji,jj+1) = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+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 1693 CALL lbc_lnk( risfdep,'T', 1. ) 1694 CALL lbc_lnk( bathy, 'T', 1. ) 1695 1696 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1697 CALL lbc_lnk( zbathy, 'T', 1. ) 1698 mbathy(:,:) = INT( zbathy(:,:) ) 1699 ENDIF 1700 ! if not compatible after all check, mask T 1701 DO jj = 1, jpj 1702 DO ji = 1, jpi 1703 IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN 1704 misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0._wp ; 1705 END IF 1706 END DO 1707 END DO 1708 1709 WHERE (mbathy(:,:) == 1) 1710 mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp 1711 END WHERE 1712 END DO 1713 ! end check compatibility ice shelf/bathy 1714 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1715 WHERE (risfdep(:,:) <= 10._wp) 1716 misfdep = 1; risfdep = 0.0_wp; 1717 END WHERE 1718 1719 IF( icompt == 0 ) THEN 1720 IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry' 1721 ELSE 1722 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry' 1723 ENDIF 1724 1725 ! compute scale factor and depth at T- and W- points 998 1726 DO jj = 1, jpj 999 1727 DO ji = 1, jpi … … 1017 1745 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1018 1746 ENDIF 1747 ! gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1019 1748 !gm Bug? check the gdepw_1d 1020 1749 ! ... on ik … … 1022 1751 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1023 1752 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1024 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 1025 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 1026 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 1027 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 1753 e3t_0 (ji,jj,ik ) = gdepw_0(ji,jj,ik+1) - gdepw_1d(ik ) 1754 e3w_0 (ji,jj,ik ) = gdept_0(ji,jj,ik ) - gdept_1d(ik-1) 1028 1755 ! ... on ik+1 1029 1756 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1030 1757 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1031 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)1032 1758 ENDIF 1033 1759 ENDIF … … 1055 1781 END DO 1056 1782 ! 1057 IF ( ln_isfcav ) THEN1058 1783 ! (ISF) Definition of e3t, u, v, w for ISF case 1059 1060 1061 1062 1063 1064 1784 DO jj = 1, jpj 1785 DO ji = 1, jpi 1786 ik = misfdep(ji,jj) 1787 IF( ik > 1 ) THEN ! ice shelf point only 1788 IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik) 1789 gdepw_0(ji,jj,ik) = risfdep(ji,jj) 1065 1790 !gm Bug? check the gdepw_0 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 e3w_0 (ji,jj,ik ) =2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))1078 1791 ! ... on ik 1792 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) & 1793 & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) & 1794 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) 1795 e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 1796 e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 1797 1798 IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column) 1799 e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 1800 ENDIF 1801 ! ... on ik / ik-1 1802 e3w_0 (ji,jj,ik ) = e3t_0 (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 1803 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 1079 1804 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 1080 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1805 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1806 ENDIF 1807 END DO 1808 END DO 1809 1810 it = 0 1811 DO jj = 1, jpj 1812 DO ji = 1, jpi 1813 ik = misfdep(ji,jj) 1814 IF( ik > 1 ) THEN ! ice shelf point only 1815 e3tp (ji,jj) = e3t_0(ji,jj,ik ) 1816 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 1817 ! test 1818 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik ) 1819 IF( zdiff <= 0. .AND. lwp ) THEN 1820 it = it + 1 1821 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1822 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 1823 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1824 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj) 1081 1825 ENDIF 1082 END DO1826 ENDIF 1083 1827 END DO 1084 !1085 it = 01086 DO jj = 1, jpj1087 DO ji = 1, jpi1088 ik = misfdep(ji,jj)1089 IF( ik > 1 ) THEN ! ice shelf point only1090 e3tp (ji,jj) = e3t_0(ji,jj,ik )1091 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )1092 ! test1093 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik )1094 IF( zdiff <= 0. .AND. lwp ) THEN1095 it = it + 11096 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj1097 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)1098 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff1099 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj)1100 ENDIF1101 ENDIF1102 END DO1103 END DO1104 END IF1105 ! END (ISF)1106 1107 ! Scale factors and depth at U-, V-, UW and VW-points1108 DO jk = 1, jpk ! initialisation to z-scale factors1109 e3u_0 (:,:,jk) = e3t_1d(jk)1110 e3v_0 (:,:,jk) = e3t_1d(jk)1111 e3uw_0(:,:,jk) = e3w_1d(jk)1112 e3vw_0(:,:,jk) = e3w_1d(jk)1113 END DO1114 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors1115 DO jj = 1, jpjm11116 DO ji = 1, fs_jpim1 ! vector opt.1117 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )1118 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )1119 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )1120 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )1121 END DO1122 END DO1123 END DO1124 IF ( ln_isfcav ) THEN1125 ! (ISF) define e3uw (adapted for 2 cells in the water column)1126 DO jj = 2, jpjm11127 DO ji = 2, fs_jpim1 ! vector opt.1128 ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj))1129 ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj))1130 IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji+1,jj ,ikb ) ) &1131 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj ,ikb-1) )1132 ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1))1133 ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1))1134 IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji ,jj+1,ikb ) ) &1135 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji ,jj+1,ikb-1) )1136 END DO1137 END DO1138 END IF1139 1140 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions1141 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp )1142 !1143 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)1144 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk)1145 WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk)1146 WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk)1147 WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk)1148 END DO1149 1150 ! Scale factor at F-point1151 DO jk = 1, jpk ! initialisation to z-scale factors1152 e3f_0(:,:,jk) = e3t_1d(jk)1153 END DO1154 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors1155 DO jj = 1, jpjm11156 DO ji = 1, fs_jpim1 ! vector opt.1157 e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )1158 END DO1159 END DO1160 END DO1161 CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions1162 !1163 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)1164 WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk)1165 END DO1166 !!gm bug ? : must be a do loop with mj0,mj11167 !1168 e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 21169 e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)1170 e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)1171 e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)1172 e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)1173 1174 ! Control of the sign1175 IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' )1176 IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' )1177 IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' )1178 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' )1179 1180 ! Compute gde3w_0 (vertical sum of e3w)1181 IF ( ln_isfcav ) THEN ! if cavity1182 WHERE( misfdep == 0 ) misfdep = 11183 DO jj = 1,jpj1184 DO ji = 1,jpi1185 gde3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)1186 DO jk = 2, misfdep(ji,jj)1187 gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)1188 END DO1189 IF( misfdep(ji,jj) >= 2 ) gde3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj))1190 DO jk = misfdep(ji,jj) + 1, jpk1191 gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)1192 END DO1193 END DO1194 END DO1195 ELSE ! no cavity1196 gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)1197 DO jk = 2, jpk1198 gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk)1199 END DO1200 END IF1201 !1202 CALL wrk_dealloc( jpi,jpj,jpk, zprt )1203 !1204 IF( nn_timing == 1 ) CALL timing_stop('zgr_zps')1205 !1206 END SUBROUTINE zgr_zps1207 1208 1209 SUBROUTINE zgr_isf1210 !!----------------------------------------------------------------------1211 !! *** ROUTINE zgr_isf ***1212 !!1213 !! ** Purpose : check the bathymetry in levels1214 !!1215 !! ** Method : THe water column have to contained at least 2 cells1216 !! Bathymetry and isfdraft are modified (dig/close) to respect1217 !! this criterion.1218 !!1219 !!1220 !! ** Action : - test compatibility between isfdraft and bathy1221 !! - bathy and isfdraft are modified1222 !!----------------------------------------------------------------------1223 INTEGER :: ji, jj, jk, jl ! dummy loop indices1224 INTEGER :: ik, it ! temporary integers1225 INTEGER :: id, jd, nprocd1226 INTEGER :: icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1 ! (ISF)1227 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points1228 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t1229 REAL(wp) :: zmax, zmin ! Maximum and minimum depth1230 REAL(wp) :: zdiff ! temporary scalar1231 REAL(wp) :: zrefdep ! temporary scalar1232 REAL(wp) :: zbathydiff, zrisfdepdiff ! isf temporary scalar1233 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH)1234 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH)1235 !!---------------------------------------------------------------------1236 !1237 IF( nn_timing == 1 ) CALL timing_start('zgr_isf')1238 !1239 CALL wrk_alloc( jpi,jpj, zbathy, zmask, zrisfdep)1240 CALL wrk_alloc( jpi,jpj, zmisfdep, zmbathy )1241 1242 1243 ! (ISF) compute misfdep1244 WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ; misfdep(:,:) = 1 ! open water : set misfdep to 11245 ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level1246 END WHERE1247 1248 ! Compute misfdep for ocean points (i.e. first wet level)1249 ! find the first ocean level such that the first level thickness1250 ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where1251 ! e3t_0 is the reference level thickness1252 DO jk = 2, jpkm11253 zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )1254 WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+11255 1828 END DO 1256 WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) > 0._wp)1257 risfdep(:,:) = 0. ; misfdep(:,:) = 11258 END WHERE1259 1260 ! 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 situation1261 icompt = 01262 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together1263 DO jl = 1, 101264 WHERE (bathy(:,:) == risfdep(:,:) )1265 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp1266 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp1267 END WHERE1268 WHERE (mbathy(:,:) <= 0)1269 misfdep(:,:) = 0; risfdep(:,:) = 0._wp1270 mbathy (:,:) = 0; bathy (:,:) = 0._wp1271 END WHERE1272 IF( lk_mpp ) THEN1273 zbathy(:,:) = FLOAT( misfdep(:,:) )1274 CALL lbc_lnk( zbathy, 'T', 1. )1275 misfdep(:,:) = INT( zbathy(:,:) )1276 CALL lbc_lnk( risfdep, 'T', 1. )1277 CALL lbc_lnk( bathy, 'T', 1. )1278 zbathy(:,:) = FLOAT( mbathy(:,:) )1279 CALL lbc_lnk( zbathy, 'T', 1. )1280 mbathy(:,:) = INT( zbathy(:,:) )1281 ENDIF1282 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN1283 misfdep( 1 ,:) = misfdep(jpim1,:) ! local domain is cyclic east-west1284 misfdep(jpi,:) = misfdep( 2 ,:)1285 ENDIF1286 1287 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN1288 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west1289 mbathy(jpi,:) = mbathy( 2 ,:)1290 ENDIF1291 1292 ! split last cell if possible (only where water column is 2 cell or less)1293 DO jk = jpkm1, 1, -11294 zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )1295 WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy)1296 mbathy(:,:) = jk1297 bathy(:,:) = zmax1298 END WHERE1299 END DO1300 1301 ! split top cell if possible (only where water column is 2 cell or less)1302 DO jk = 2, jpkm11303 zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )1304 WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy)1305 misfdep(:,:) = jk1306 risfdep(:,:) = zmax1307 END WHERE1308 END DO1309 1310 1311 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition1312 DO jj = 1, jpj1313 DO ji = 1, jpi1314 ! find the minimum change option:1315 ! test bathy1316 IF (risfdep(ji,jj) > 1) THEN1317 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) &1318 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))1319 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) &1320 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))1321 1322 IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN1323 IF (zbathydiff .LE. zrisfdepdiff) THEN1324 bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )1325 mbathy(ji,jj)= mbathy(ji,jj) + 11326 ELSE1327 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )1328 misfdep(ji,jj) = misfdep(ji,jj) - 11329 END IF1330 END IF1331 END IF1332 END DO1333 END DO1334 1335 ! At least 2 levels for water thickness at T, U, and V point.1336 DO jj = 1, jpj1337 DO ji = 1, jpi1338 ! find the minimum change option:1339 ! test bathy1340 IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN1341 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1)&1342 & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))1343 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) &1344 & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))1345 IF (zbathydiff .LE. zrisfdepdiff) THEN1346 mbathy(ji,jj) = mbathy(ji,jj) + 11347 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )1348 ELSE1349 misfdep(ji,jj)= misfdep(ji,jj) - 11350 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )1351 END IF1352 ENDIF1353 END DO1354 END DO1355 1356 ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1)1357 DO jj = 1, jpjm11358 DO ji = 1, jpim11359 IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN1360 zbathydiff =ABS(bathy(ji,jj ) - (gdepw_1d(mbathy (ji,jj)+1) &1361 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat )))1362 zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1)) &1363 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat )))1364 IF (zbathydiff .LE. zrisfdepdiff) THEN1365 mbathy(ji,jj) = mbathy(ji,jj) + 11366 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) &1367 & + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat )1368 ELSE1369 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 11370 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) &1371 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )1372 END IF1373 ENDIF1374 END DO1375 END DO1376 1377 IF( lk_mpp ) THEN1378 zbathy(:,:) = FLOAT( misfdep(:,:) )1379 CALL lbc_lnk( zbathy, 'T', 1. )1380 misfdep(:,:) = INT( zbathy(:,:) )1381 CALL lbc_lnk( risfdep, 'T', 1. )1382 CALL lbc_lnk( bathy, 'T', 1. )1383 zbathy(:,:) = FLOAT( mbathy(:,:) )1384 CALL lbc_lnk( zbathy, 'T', 1. )1385 mbathy(:,:) = INT( zbathy(:,:) )1386 ENDIF1387 ! point V misdep(ji,jj) EQ mbathy(ji,jj+1)1388 DO jj = 1, jpjm11389 DO ji = 1, jpim11390 IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN1391 zbathydiff =ABS( bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) &1392 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )))1393 zrisfdepdiff=ABS(risfdep(ji,jj ) - (gdepw_1d(misfdep(ji,jj ) ) &1394 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat )))1395 IF (zbathydiff .LE. zrisfdepdiff) THEN1396 mbathy (ji,jj+1) = mbathy(ji,jj+1) + 11397 bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )1398 ELSE1399 misfdep(ji,jj) = misfdep(ji,jj) - 11400 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat )1401 END IF1402 ENDIF1403 END DO1404 END DO1405 1406 1407 IF( lk_mpp ) THEN1408 zbathy(:,:) = FLOAT( misfdep(:,:) )1409 CALL lbc_lnk( zbathy, 'T', 1. )1410 misfdep(:,:) = INT( zbathy(:,:) )1411 CALL lbc_lnk( risfdep, 'T', 1. )1412 CALL lbc_lnk( bathy, 'T', 1. )1413 zbathy(:,:) = FLOAT( mbathy(:,:) )1414 CALL lbc_lnk( zbathy, 'T', 1. )1415 mbathy(:,:) = INT( zbathy(:,:) )1416 ENDIF1417 1418 ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1)1419 DO jj = 1, jpjm11420 DO ji = 1, jpim11421 IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN1422 zbathydiff =ABS( bathy(ji ,jj) - (gdepw_1d(mbathy (ji,jj)+1) &1423 & + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat )))1424 zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) &1425 & - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat )))1426 IF (zbathydiff .LE. zrisfdepdiff) THEN1427 mbathy(ji,jj) = mbathy(ji,jj) + 11428 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )1429 ELSE1430 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 11431 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )1432 END IF1433 ENDIF1434 ENDDO1435 ENDDO1436 1437 IF( lk_mpp ) THEN1438 zbathy(:,:) = FLOAT( misfdep(:,:) )1439 CALL lbc_lnk( zbathy, 'T', 1. )1440 misfdep(:,:) = INT( zbathy(:,:) )1441 CALL lbc_lnk( risfdep, 'T', 1. )1442 CALL lbc_lnk( bathy, 'T', 1. )1443 zbathy(:,:) = FLOAT( mbathy(:,:) )1444 CALL lbc_lnk( zbathy, 'T', 1. )1445 mbathy(:,:) = INT( zbathy(:,:) )1446 ENDIF1447 1448 ! point U misfdep(ji,jj) EQ bathy(ji,jj+1)1449 DO jj = 1, jpjm11450 DO ji = 1, jpim11451 IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN1452 zbathydiff =ABS( bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) &1453 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat )))1454 zrisfdepdiff=ABS(risfdep(ji ,jj) - (gdepw_1d(misfdep(ji ,jj) ) &1455 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat )))1456 IF (zbathydiff .LE. zrisfdepdiff) THEN1457 mbathy(ji+1,jj) = mbathy (ji+1,jj) + 11458 bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) &1459 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat )1460 ELSE1461 misfdep(ji,jj) = misfdep(ji ,jj) - 11462 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) &1463 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat )1464 END IF1465 ENDIF1466 ENDDO1467 ENDDO1468 1469 IF( lk_mpp ) THEN1470 zbathy(:,:) = FLOAT( misfdep(:,:) )1471 CALL lbc_lnk( zbathy, 'T', 1. )1472 misfdep(:,:) = INT( zbathy(:,:) )1473 CALL lbc_lnk( risfdep, 'T', 1. )1474 CALL lbc_lnk( bathy, 'T', 1. )1475 zbathy(:,:) = FLOAT( mbathy(:,:) )1476 CALL lbc_lnk( zbathy, 'T', 1. )1477 mbathy(:,:) = INT( zbathy(:,:) )1478 ENDIF1479 END DO1480 ! end dig bathy/ice shelf to be compatible1481 ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness1482 DO jl = 1,201483 1484 ! remove single point "bay" on isf coast line in the ice shelf draft'1485 DO jk = 2, jpk1486 WHERE (misfdep==0) misfdep=jpk1487 zmask=01488 WHERE (misfdep .LE. jk) zmask=11489 DO jj = 2, jpjm11490 DO ji = 2, jpim11491 IF (misfdep(ji,jj) .EQ. jk) THEN1492 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)1493 IF (ibtest .LE. 1) THEN1494 risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+11495 IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) misfdep(ji,jj) = jpk1496 END IF1497 END IF1498 END DO1499 END DO1500 END DO1501 WHERE (misfdep==jpk)1502 misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.1503 END WHERE1504 IF( lk_mpp ) THEN1505 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 ENDIF1514 1515 ! remove single point "bay" on bathy coast line beneath an ice shelf'1516 DO jk = jpk,1,-11517 zmask=01518 WHERE (mbathy .GE. jk ) zmask=11519 DO jj = 2, jpjm11520 DO ji = 2, jpim11521 IF (mbathy(ji,jj) .EQ. jk .AND. misfdep(ji,jj) .GE. 2) THEN1522 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)1523 IF (ibtest .LE. 1) THEN1524 bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-11525 IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) mbathy(ji,jj) = 01526 END IF1527 END IF1528 END DO1529 END DO1530 END DO1531 WHERE (mbathy==0)1532 misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.1533 END WHERE1534 IF( lk_mpp ) THEN1535 zbathy(:,:) = FLOAT( misfdep(:,:) )1536 CALL lbc_lnk( zbathy, 'T', 1. )1537 misfdep(:,:) = INT( zbathy(:,:) )1538 CALL lbc_lnk( risfdep, 'T', 1. )1539 CALL lbc_lnk( bathy, 'T', 1. )1540 zbathy(:,:) = FLOAT( mbathy(:,:) )1541 CALL lbc_lnk( zbathy, 'T', 1. )1542 mbathy(:,:) = INT( zbathy(:,:) )1543 ENDIF1544 1545 ! fill hole in ice shelf1546 zmisfdep = misfdep1547 zrisfdep = risfdep1548 WHERE (zmisfdep .LE. 1) zmisfdep=jpk1549 DO jj = 2, jpjm11550 DO ji = 2, jpim11551 ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj )1552 ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1)1553 IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj ) - 1)1554 IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj ) - 1)1555 IF( zmisfdep(ji,jj) .GE. mbathy(ji ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji ,jj-1) - 1)1556 IF( zmisfdep(ji,jj) .GE. mbathy(ji ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji ,jj+1) - 1)1557 ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)1558 IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN1559 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp1560 END IF1561 IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) .GE. 2) THEN1562 misfdep(ji,jj) = ibtest1563 risfdep(ji,jj) = gdepw_1d(ibtest)1564 ENDIF1565 ENDDO1566 ENDDO1567 1568 IF( lk_mpp ) THEN1569 zbathy(:,:) = FLOAT( misfdep(:,:) )1570 CALL lbc_lnk( zbathy, 'T', 1. )1571 misfdep(:,:) = INT( zbathy(:,:) )1572 CALL lbc_lnk( risfdep, 'T', 1. )1573 CALL lbc_lnk( bathy, 'T', 1. )1574 zbathy(:,:) = FLOAT( mbathy(:,:) )1575 CALL lbc_lnk( zbathy, 'T', 1. )1576 mbathy(:,:) = INT( zbathy(:,:) )1577 ENDIF1578 !1579 !! fill hole in bathymetry1580 zmbathy (:,:)=mbathy (:,:)1581 DO jj = 2, jpjm11582 DO ji = 2, jpim11583 ibtestim1 = zmbathy(ji-1,jj ) ; ibtestip1 = zmbathy(ji+1,jj )1584 ibtestjm1 = zmbathy(ji ,jj-1) ; ibtestjp1 = zmbathy(ji ,jj+1)1585 IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj ) + 1)1586 IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj ) ) ibtestip1 = 01587 IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj-1) ) ibtestjm1 = 01588 IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj+1) ) ibtestjp1 = 01589 ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)1590 IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN1591 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ;1592 END IF1593 IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) .GE. 2) THEN1594 mbathy(ji,jj) = ibtest1595 bathy(ji,jj) = gdepw_1d(ibtest+1)1596 ENDIF1597 END DO1598 END DO1599 IF( lk_mpp ) THEN1600 zbathy(:,:) = FLOAT( misfdep(:,:) )1601 CALL lbc_lnk( zbathy, 'T', 1. )1602 misfdep(:,:) = INT( zbathy(:,:) )1603 CALL lbc_lnk( risfdep, 'T', 1. )1604 CALL lbc_lnk( bathy, 'T', 1. )1605 zbathy(:,:) = FLOAT( mbathy(:,:) )1606 CALL lbc_lnk( zbathy, 'T', 1. )1607 mbathy(:,:) = INT( zbathy(:,:) )1608 ENDIF1609 ! if not compatible after all check (ie U point water column less than 2 cells), mask U1610 DO jj = 1, jpjm11611 DO ji = 1, jpim11612 IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN1613 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ;1614 END IF1615 END DO1616 END DO1617 IF( lk_mpp ) THEN1618 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 ENDIF1627 ! if not compatible after all check (ie U point water column less than 2 cells), mask U1628 DO jj = 1, jpjm11629 DO ji = 1, jpim11630 IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN1631 mbathy(ji+1,jj) = mbathy(ji+1,jj) - 1; bathy(ji+1,jj) = gdepw_1d(mbathy(ji+1,jj)+1) ;1632 END IF1633 END DO1634 END DO1635 IF( lk_mpp ) THEN1636 zbathy(:,:) = FLOAT( misfdep(:,:) )1637 CALL lbc_lnk( zbathy, 'T', 1. )1638 misfdep(:,:) = INT( zbathy(:,:) )1639 CALL lbc_lnk( risfdep, 'T', 1. )1640 CALL lbc_lnk( bathy, 'T', 1. )1641 zbathy(:,:) = FLOAT( mbathy(:,:) )1642 CALL lbc_lnk( zbathy, 'T', 1. )1643 mbathy(:,:) = INT( zbathy(:,:) )1644 ENDIF1645 ! if not compatible after all check (ie V point water column less than 2 cells), mask V1646 DO jj = 1, jpjm11647 DO ji = 1, jpi1648 IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN1649 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ;1650 END IF1651 END DO1652 END DO1653 IF( lk_mpp ) THEN1654 zbathy(:,:) = FLOAT( misfdep(:,:) )1655 CALL lbc_lnk( zbathy, 'T', 1. )1656 misfdep(:,:) = INT( zbathy(:,:) )1657 CALL lbc_lnk( risfdep, 'T', 1. )1658 CALL lbc_lnk( bathy, 'T', 1. )1659 zbathy(:,:) = FLOAT( mbathy(:,:) )1660 CALL lbc_lnk( zbathy, 'T', 1. )1661 mbathy(:,:) = INT( zbathy(:,:) )1662 ENDIF1663 ! if not compatible after all check (ie V point water column less than 2 cells), mask V1664 DO jj = 1, jpjm11665 DO ji = 1, jpi1666 IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN1667 mbathy(ji,jj+1) = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ;1668 END IF1669 END DO1670 END DO1671 IF( lk_mpp ) THEN1672 zbathy(:,:) = FLOAT( misfdep(:,:) )1673 CALL lbc_lnk( zbathy, 'T', 1. )1674 misfdep(:,:) = INT( zbathy(:,:) )1675 CALL lbc_lnk( risfdep, 'T', 1. )1676 CALL lbc_lnk( bathy, 'T', 1. )1677 zbathy(:,:) = FLOAT( mbathy(:,:) )1678 CALL lbc_lnk( zbathy, 'T', 1. )1679 mbathy(:,:) = INT( zbathy(:,:) )1680 ENDIF1681 ! if not compatible after all check, mask T1682 DO jj = 1, jpj1683 DO ji = 1, jpi1684 IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN1685 misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0._wp ;1686 END IF1687 END DO1688 END DO1689 1690 WHERE (mbathy(:,:) == 1)1691 mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp1692 END WHERE1693 END DO1694 ! end check compatibility ice shelf/bathy1695 ! remove very shallow ice shelf (less than ~ 10m if 75L)1696 WHERE (misfdep(:,:) <= 5)1697 misfdep = 1; risfdep = 0.0_wp;1698 END WHERE1699 1700 IF( icompt == 0 ) THEN1701 IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry'1702 ELSE1703 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry'1704 ENDIF1705 1829 1706 1830 CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) … … 1709 1833 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf') 1710 1834 ! 1711 END SUBROUTINE 1835 END SUBROUTINE zgr_isf 1712 1836 1713 1837
Note: See TracChangeset
for help on using the changeset viewer.
