- Timestamp:
- 2015-07-20T19:43:15+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r5506 r5619 365 365 !! - bathy : meter bathymetry (in meters) 366 366 !!---------------------------------------------------------------------- 367 INTEGER :: ji, jj, j l, jk ! dummy loop indices367 INTEGER :: ji, jj, jk ! dummy loop indices 368 368 INTEGER :: inum ! temporary logical unit 369 369 INTEGER :: ierror ! error flag … … 472 472 risfdep(:,:)=0.e0 473 473 misfdep(:,:)=1 474 ! 475 ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code 476 IF( cp_cfg == "isomip" .AND. ln_isfcav ) THEN 477 risfdep(:,:)=200.e0 478 misfdep(:,:)=1 479 ij0 = 1 ; ij1 = 40 480 DO jj = mj0(ij0), mj1(ij1) 481 risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 482 END DO 483 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 484 ! 485 ELSEIF ( cp_cfg == "isomip2" .AND. ln_isfcav ) THEN 486 ! 487 risfdep(:,:)=0.e0 488 misfdep(:,:)=1 489 ij0 = 1 ; ij1 = 40 490 DO jj = mj0(ij0), mj1(ij1) 491 risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 492 END DO 493 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 494 END IF 474 495 ! 475 496 DEALLOCATE( idta, zdta ) … … 529 550 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 530 551 END IF 552 ! set grounded point to 0 553 WHERE (bathy(:,:) .LE. risfdep(:,:)+1e-2 ) 554 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 555 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 556 END WHERE 531 557 ! 532 558 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration … … 952 978 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 953 979 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t 954 REAL(wp) :: zmax ! Maximum depth955 980 REAL(wp) :: zdiff ! temporary scalar 956 REAL(wp) :: z refdep! temporary scalar981 REAL(wp) :: zmax ! temporary scalar 957 982 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt 958 983 !!--------------------------------------------------------------------- … … 993 1018 END DO 994 1019 995 IF ( ln_isfcav ) CALL zgr_isf996 997 1020 ! Scale factors and depth at T- and W-points 998 1021 DO jk = 1, jpk ! intitialization to the reference z-coordinate … … 1002 1025 e3w_0 (:,:,jk) = e3w_1d (jk) 1003 1026 END DO 1004 ! 1005 DO jj = 1, jpj 1006 DO ji = 1, jpi 1007 ik = mbathy(ji,jj) 1008 IF( ik > 0 ) THEN ! ocean point only 1009 ! max ocean level case 1010 IF( ik == jpkm1 ) THEN 1011 zdepwp = bathy(ji,jj) 1012 ze3tp = bathy(ji,jj) - gdepw_1d(ik) 1013 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 1014 e3t_0(ji,jj,ik ) = ze3tp 1015 e3t_0(ji,jj,ik+1) = ze3tp 1016 e3w_0(ji,jj,ik ) = ze3wp 1017 e3w_0(ji,jj,ik+1) = ze3tp 1018 gdepw_0(ji,jj,ik+1) = zdepwp 1019 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 1020 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 1021 ! 1022 ELSE ! standard case 1023 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 1024 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1027 1028 ! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf 1029 IF ( ln_isfcav ) CALL zgr_isf 1030 1031 ! Scale factors and depth at T- and W-points 1032 IF ( .NOT. ln_isfcav ) THEN 1033 DO jj = 1, jpj 1034 DO ji = 1, jpi 1035 ik = mbathy(ji,jj) 1036 IF( ik > 0 ) THEN ! ocean point only 1037 ! max ocean level case 1038 IF( ik == jpkm1 ) THEN 1039 zdepwp = bathy(ji,jj) 1040 ze3tp = bathy(ji,jj) - gdepw_1d(ik) 1041 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 1042 e3t_0(ji,jj,ik ) = ze3tp 1043 e3t_0(ji,jj,ik+1) = ze3tp 1044 e3w_0(ji,jj,ik ) = ze3wp 1045 e3w_0(ji,jj,ik+1) = ze3tp 1046 gdepw_0(ji,jj,ik+1) = zdepwp 1047 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 1048 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 1049 ! 1050 ELSE ! standard case 1051 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 1052 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1053 ENDIF 1054 !gm Bug? check the gdepw_1d 1055 ! ... on ik 1056 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 1057 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1058 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1059 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 1060 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 1061 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 1062 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 1063 ! ... on ik+1 1064 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1065 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1066 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 1025 1067 ENDIF 1026 !gm Bug? check the gdepw_1d1027 ! ... on ik1028 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) )1033 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) &1034 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )1035 ! ... on ik+11036 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)1037 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)1038 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)1039 1068 ENDIF 1040 END IF1041 END DO 1042 END DO1043 !1044 it = 01045 DO jj = 1, jpj1046 DO ji = 1, jpi1047 ik = mbathy(ji,jj)1048 IF( ik > 0 ) THEN ! ocean point only1049 e3tp (ji,jj) = e3t_0(ji,jj,ik)1050 e3wp (ji,jj) = e3w_0(ji,jj,ik)1051 ! test1052 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik )1053 IF( zdiff <= 0._wp .AND. lwp ) THEN1054 it = it + 11055 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj1056 WRITE(numout,*) ' bathy = ', bathy(ji,jj)1057 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff1058 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik )1069 END DO 1070 END DO 1071 ! 1072 it = 0 1073 DO jj = 1, jpj 1074 DO ji = 1, jpi 1075 ik = mbathy(ji,jj) 1076 IF( ik > 0 ) THEN ! ocean point only 1077 e3tp (ji,jj) = e3t_0(ji,jj,ik) 1078 e3wp (ji,jj) = e3w_0(ji,jj,ik) 1079 ! test 1080 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik ) 1081 IF( zdiff <= 0._wp .AND. lwp ) THEN 1082 it = it + 1 1083 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1084 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 1085 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1086 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik ) 1087 ENDIF 1059 1088 ENDIF 1060 ENDIF 1061 END DO 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 1089 END DO 1090 END DO 1111 1091 END IF 1112 ! END (ISF) 1113 1092 ! 1114 1093 ! Scale factors and depth at U-, V-, UW and VW-points 1115 1094 DO jk = 1, jpk ! initialisation to z-scale factors … … 1119 1098 e3vw_0(:,:,jk) = e3w_1d(jk) 1120 1099 END DO 1100 1121 1101 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors 1122 1102 DO jj = 1, jpjm1 … … 1148 1128 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 1149 1129 ! 1130 1150 1131 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1151 1132 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk) … … 1255 1236 !!---------------------------------------------------------------------- 1256 1237 !! 1257 INTEGER :: ji, jj, jk, jl ! dummy loop indices 1258 INTEGER :: ik, it ! temporary integers 1259 INTEGER :: id, jd, nprocd 1238 INTEGER :: ji, jj, jl, jk ! dummy loop indices 1239 INTEGER :: ik, it ! temporary integers 1260 1240 INTEGER :: icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1 ! (ISF) 1261 LOGICAL :: ll_print ! Allow control print for debugging 1241 REAL(wp) :: zdepth ! Ajusted ocean depth to avoid too small e3t 1242 REAL(wp) :: zmax ! Maximum and minimum depth 1243 REAL(wp) :: zbathydiff, zrisfdepdiff ! isf temporary scalar 1262 1244 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 1245 REAL(wp) :: zdepwp ! Ajusted ocean depth to avoid too small e3t 1265 1246 REAL(wp) :: zdiff ! temporary scalar 1266 REAL(wp) :: zrefdep ! temporary scalar1267 REAL(wp) :: zbathydiff, zrisfdepdiff ! isf temporary scalar1268 1247 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH) 1269 1248 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH) … … 1280 1259 ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level 1281 1260 END WHERE 1261 1262 ! set grounded point to 0 1263 WHERE (bathy(:,:) .LE. risfdep(:,:)+1e-2 ) 1264 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 1265 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 1266 END WHERE 1282 1267 1283 1268 ! Compute misfdep for ocean points (i.e. first wet level) … … 1292 1277 risfdep(:,:) = 0. ; misfdep(:,:) = 1 1293 1278 END WHERE 1279 1280 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1281 IF ( cp_cfg .NE. "isomip" ) THEN 1282 WHERE (risfdep(:,:) < 100 ) 1283 misfdep = 1; risfdep = 0.0_wp; 1284 END WHERE 1285 END IF 1294 1286 1295 1287 ! 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 … … 1297 1289 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 1298 1290 DO jl = 1, 10 1299 WHERE (bathy(:,:) . EQ. risfdep(:,:))1291 WHERE (bathy(:,:) .LE. risfdep(:,:)+1e-2 ) 1300 1292 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 1301 1293 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp … … 1326 1318 1327 1319 ! split last cell if possible (only where water column is 2 cell or less) 1328 DO jk = jpkm1, 1, -11329 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(:,:) = jk1332 bathy(:,:) = zmax1333 END WHERE1334 END DO1320 !DO jk = jpkm1, 1, -1 1321 ! zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1322 ! WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 1323 ! mbathy(:,:) = jk 1324 ! bathy(:,:) = zmax 1325 ! END WHERE 1326 !END DO 1335 1327 1336 1328 ! split top cell if possible (only where water column is 2 cell or less) … … 1350 1342 ! test bathy 1351 1343 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 ))) 1344 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1345 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1356 1346 1357 1347 IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN 1358 IF (zbathydiff .LE. zrisfdepdiff) THEN1359 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) + 11361 ELSE1348 ! IF (zbathydiff .LE. zrisfdepdiff) THEN 1349 ! bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 1350 ! mbathy(ji,jj)= mbathy(ji,jj) + 1 1351 ! ELSE 1362 1352 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 1363 1353 misfdep(ji,jj) = misfdep(ji,jj) - 1 1364 END IF1354 ! END IF 1365 1355 END IF 1366 1356 END IF … … 1374 1364 ! test bathy 1375 1365 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 1366 ! zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1367 ! zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1368 ! IF (zbathydiff .LE. zrisfdepdiff) THEN 1369 ! mbathy(ji,jj) = mbathy(ji,jj) + 1 1370 ! bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1371 ! ELSE 1384 1372 misfdep(ji,jj)= misfdep(ji,jj) - 1 1385 1373 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 1386 END IF1374 ! END IF 1387 1375 ENDIF 1388 1376 END DO … … 1393 1381 DO ji = 1, jpim1 1394 1382 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 1383 ! zbathydiff =ABS(bathy(ji,jj ) - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat ))) 1384 ! zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 1385 ! IF (zbathydiff .LE. zrisfdepdiff) THEN 1386 ! mbathy(ji,jj) = mbathy(ji,jj) + 1 1387 ! bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat ) 1388 ! ELSE 1404 1389 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 1390 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 1391 ! END IF 1408 1392 ENDIF 1409 1393 END DO … … 1424 1408 DO ji = 1, jpim1 1425 1409 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 1410 ! zbathydiff =ABS( bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 1411 ! zrisfdepdiff=ABS(risfdep(ji,jj ) - (gdepw_1d(misfdep(ji,jj ) ) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat ))) 1412 ! IF (zbathydiff .LE. zrisfdepdiff) THEN 1413 ! mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 1414 ! bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 1415 ! ELSE 1434 1416 misfdep(ji,jj) = misfdep(ji,jj) - 1 1435 1417 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat ) 1436 END IF1418 ! END IF 1437 1419 ENDIF 1438 1420 END DO … … 1455 1437 DO ji = 1, jpim1 1456 1438 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 1439 ! zbathydiff =ABS( bathy(ji ,jj) - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat ))) 1440 ! zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 1441 ! IF (zbathydiff .LE. zrisfdepdiff) THEN 1442 ! mbathy(ji,jj) = mbathy(ji,jj) + 1 1443 ! bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1444 ! ELSE 1465 1445 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 1466 1446 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 1467 END IF1447 ! END IF 1468 1448 ENDIF 1469 1449 ENDDO … … 1485 1465 DO ji = 1, jpim1 1486 1466 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 1467 ! zbathydiff =ABS( bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 1468 ! zrisfdepdiff=ABS(risfdep(ji ,jj) - (gdepw_1d(misfdep(ji ,jj) ) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat ))) 1469 ! IF (zbathydiff .LE. zrisfdepdiff) THEN 1470 ! mbathy(ji+1,jj) = mbathy (ji+1,jj) + 1 1471 ! bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 1472 ! ELSE 1496 1473 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 1474 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat ) 1475 ! END IF 1500 1476 ENDIF 1501 1477 ENDDO … … 1729 1705 ! end check compatibility ice shelf/bathy 1730 1706 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1731 WHERE (misfdep(:,:) <= 5)1732 misfdep = 1; risfdep = 0.0_wp;1733 END WHERE1734 1735 1707 IF( icompt == 0 ) THEN 1736 1708 IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry' … … 1738 1710 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry' 1739 1711 ENDIF 1712 1713 ! compute scale factor and depth at T- and W- points 1714 DO jj = 1, jpj 1715 DO ji = 1, jpi 1716 ik = mbathy(ji,jj) 1717 IF( ik > 0 ) THEN ! ocean point only 1718 ! max ocean level case 1719 IF( ik == jpkm1 ) THEN 1720 zdepwp = bathy(ji,jj) 1721 ze3tp = bathy(ji,jj) - gdepw_1d(ik) 1722 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 1723 e3t_0(ji,jj,ik ) = ze3tp 1724 e3t_0(ji,jj,ik+1) = ze3tp 1725 e3w_0(ji,jj,ik ) = ze3wp 1726 e3w_0(ji,jj,ik+1) = ze3tp 1727 gdepw_0(ji,jj,ik+1) = zdepwp 1728 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 1729 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 1730 ! 1731 ELSE ! standard case 1732 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 1733 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1734 ENDIF 1735 ! gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1736 !gm Bug? check the gdepw_1d 1737 ! ... on ik 1738 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 1739 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1740 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1741 e3t_0 (ji,jj,ik ) = gdepw_0(ji,jj,ik+1) - gdepw_1d(ik ) 1742 e3w_0 (ji,jj,ik ) = gdept_0(ji,jj,ik ) - gdept_1d(ik-1) 1743 ! ... on ik+1 1744 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1745 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1746 ENDIF 1747 ENDIF 1748 END DO 1749 END DO 1750 ! 1751 it = 0 1752 DO jj = 1, jpj 1753 DO ji = 1, jpi 1754 ik = mbathy(ji,jj) 1755 IF( ik > 0 ) THEN ! ocean point only 1756 e3tp (ji,jj) = e3t_0(ji,jj,ik) 1757 e3wp (ji,jj) = e3w_0(ji,jj,ik) 1758 ! test 1759 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik ) 1760 IF( zdiff <= 0._wp .AND. lwp ) THEN 1761 it = it + 1 1762 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1763 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 1764 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1765 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik ) 1766 ENDIF 1767 ENDIF 1768 END DO 1769 END DO 1770 ! 1771 ! (ISF) Definition of e3t, u, v, w for ISF case 1772 DO jj = 1, jpj 1773 DO ji = 1, jpi 1774 ik = misfdep(ji,jj) 1775 IF( ik > 1 ) THEN ! ice shelf point only 1776 IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik) 1777 gdepw_0(ji,jj,ik) = risfdep(ji,jj) 1778 !gm Bug? check the gdepw_0 1779 ! ... on ik 1780 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) & 1781 & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) & 1782 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) 1783 e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 1784 e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 1785 1786 IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column) 1787 e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 1788 ENDIF 1789 ! ... on ik / ik-1 1790 e3w_0 (ji,jj,ik ) = e3t_0 (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 1791 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 1792 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 1793 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1794 ENDIF 1795 END DO 1796 END DO 1797 1798 it = 0 1799 DO jj = 1, jpj 1800 DO ji = 1, jpi 1801 ik = misfdep(ji,jj) 1802 IF( ik > 1 ) THEN ! ice shelf point only 1803 e3tp (ji,jj) = e3t_0(ji,jj,ik ) 1804 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 1805 ! test 1806 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik ) 1807 IF( zdiff <= 0. .AND. lwp ) THEN 1808 it = it + 1 1809 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1810 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 1811 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1812 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj) 1813 ENDIF 1814 ENDIF 1815 END DO 1816 END DO 1740 1817 1741 1818 CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep )
Note: See TracChangeset
for help on using the changeset viewer.
