Changeset 5944 for branches/2015/dev_r5151_UKMO_ISF/NEMOGCM
- Timestamp:
- 2015-11-29T20:28:41+01:00 (8 years ago)
- Location:
- branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 13 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90
r5621 r5944 137 137 DO jk=1,nb_ana 138 138 DO ji=1,jpmax_harmo 139 IF (TRIM(tname(jk)) .eq.Wave(ji)%cname_tide) THEN139 IF (TRIM(tname(jk)) == Wave(ji)%cname_tide) THEN 140 140 name(jk) = ji 141 141 EXIT … … 362 362 X1=ana_amp(ji,jj,jh,1) 363 363 X2=-ana_amp(ji,jj,jh,2) 364 out_v(ji,jj, jh)=X1 * vmask_i(ji,jj)365 out_v(ji,jj,nb_ana+jh)=X2 * vmask_i(ji,jj)364 out_v(ji,jj, jh)=X1 * ssvmask(ji,jj) 365 out_v(ji,jj,nb_ana+jh)=X2 * ssvmask(ji,jj) 366 366 END DO 367 367 END DO … … 492 492 DO jj_sd = ji_sd, ninco 493 493 zval2 = ABS(ztmp3(ji_sd,jj_sd)) 494 IF( zval2 .GE.zval1 )THEN494 IF( zval2 >= zval1 )THEN 495 495 ipivot(ji_sd) = jj_sd 496 496 zval1 = zval2 -
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
r5624 r5944 211 211 CALL iom_put( 'bgsaline' , zdiff_sc1 / zvol_tot) ! Salt content variation (psu) 212 212 CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rau0 * rcp ) ! Heat content variation (1.e20 J) 213 CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9 )! Salt content variation (psu*km3)214 CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9) ! volume ssh variation (km3)213 CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9 ) ! Salt content variation (psu*km3) 214 CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9 ) ! volume ssh variation (km3) 215 215 CALL iom_put( 'bgfrcvol' , frc_v * 1.e-9 ) ! vol - surface forcing (km3) 216 216 CALL iom_put( 'bgfrctem' , frc_t / zvol_tot ) ! hc - surface forcing (C) -
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
r5621 r5944 319 319 WRITE(numout,*) ' cross land advection nn_cla = ', nn_cla 320 320 ENDIF 321 IF ( nn_cla .EQ.1 ) THEN321 IF ( nn_cla == 1 ) THEN 322 322 IF ( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 323 323 CONTINUE -
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r5624 r5944 332 332 ! -------------------------------------------------- 333 333 IF( ln_vvl_ztilde ) THEN 334 IF( kt .GT.nit000 ) THEN334 IF( kt > nit000 ) THEN 335 335 DO jk = 1, jpkm1 336 336 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:) & … … 426 426 IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain 427 427 ! - ML - test: for the moment, stop simulation for too large e3_t variations 428 IF( ( z_tmax .GT. rn_zdef_max ) .OR. ( z_tmin .LT.- rn_zdef_max ) ) THEN428 IF( ( z_tmax > rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN 429 429 IF( lk_mpp ) THEN 430 430 CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) ) -
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r5621 r5944 824 824 END SUBROUTINE zgr_bot_level 825 825 826 827 !!---------------------------------------------------------------------- 828 !! *** ROUTINE zgr_ bot_level ***826 SUBROUTINE zgr_top_level 827 !!---------------------------------------------------------------------- 828 !! *** ROUTINE zgr_top_level *** 829 829 !! 830 830 !! ** Purpose : defines the vertical index of ocean top (mik. arrays) … … 1149 1149 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1150 1150 END DO 1151 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))1151 IF (misfdep(ji,jj) >= 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 1152 1152 DO jk = misfdep(ji,jj) + 1, jpk 1153 1153 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) … … 1212 1212 INTEGER :: ji, jj, jl, jk ! dummy loop indices 1213 1213 INTEGER :: ik, it ! temporary integers 1214 INTEGER :: icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1 ! (ISF) 1214 INTEGER :: icompt, ibtest ! (ISF) 1215 INTEGER :: ibtestim1, ibtestip1 ! (ISF) 1216 INTEGER :: ibtestjm1, ibtestjp1 ! (ISF) 1215 1217 REAL(wp) :: zdepth ! Ajusted ocean depth to avoid too small e3t 1216 1218 REAL(wp) :: zmax ! Maximum and minimum depth 1217 REAL(wp) :: zbathydiff, zrisfdepdiff ! isf temporary scalar 1219 REAL(wp) :: zbathydiff ! isf temporary scalar 1220 REAL(wp) :: zrisfdepdiff ! isf temporary scalar 1218 1221 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 1219 1222 REAL(wp) :: zdepwp ! Ajusted ocean depth to avoid too small e3t … … 1230 1233 1231 1234 ! (ISF) compute misfdep 1232 WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ; misfdep(:,:) = 1 ! open water : set misfdep to 11233 ELSEWHERE ; 1235 WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) /= 0 ) ; misfdep(:,:) = 1 ! open water : set misfdep to 1 1236 ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level 1234 1237 END WHERE 1235 1238 … … 1242 1245 WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+1 1243 1246 END DO 1244 WHERE ( risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp)1247 WHERE ( 0._wp < risfdep(:,:) .AND. risfdep(:,:) <= e3t_1d(1) ) 1245 1248 risfdep(:,:) = 0. ; misfdep(:,:) = 1 1246 1249 END WHERE 1247 1250 1248 1251 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1249 IF ( cp_cfg .NE. "isomip" ) THEN 1250 WHERE (misfdep(:,:) <= 10 .AND. misfdep(:,:) .GT. 1) 1251 misfdep = 0; risfdep = 0.0_wp; 1252 mbathy = 0; bathy = 0.0_wp; 1253 END WHERE 1254 WHERE (bathy(:,:) <= 30.0 .AND. gphit < -60) 1255 misfdep = 0; risfdep = 0.0_wp; 1256 mbathy = 0; bathy = 0.0_wp; 1257 END WHERE 1258 END IF 1252 WHERE (risfdep(:,:) <= 10._wp .AND. misfdep(:,:) > 1) 1253 misfdep = 0; risfdep = 0.0_wp; 1254 mbathy = 0; bathy = 0.0_wp; 1255 END WHERE 1256 WHERE (bathy(:,:) <= 30.0_wp .AND. gphit < -60._wp) 1257 misfdep = 0; risfdep = 0.0_wp; 1258 mbathy = 0; bathy = 0.0_wp; 1259 END WHERE 1259 1260 1260 1261 ! 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 … … 1262 1263 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 1263 1264 DO jl = 1, 10 1264 WHERE (bathy(:,:) .EQ.risfdep(:,:) )1265 WHERE (bathy(:,:) == risfdep(:,:) ) 1265 1266 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 1266 1267 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp … … 1271 1272 ENDWHERE 1272 1273 IF( lk_mpp ) THEN 1273 zbathy(:,:) = FLOAT( misfdep(:,:) )1274 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1274 1275 CALL lbc_lnk( zbathy, 'T', 1. ) 1275 1276 misfdep(:,:) = INT( zbathy(:,:) ) 1276 CALL lbc_lnk( risfdep, 'T', 1. ) 1277 CALL lbc_lnk( bathy, 'T', 1. ) 1278 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1277 1278 CALL lbc_lnk( risfdep,'T', 1. ) 1279 CALL lbc_lnk( bathy, 'T', 1. ) 1280 1281 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1279 1282 CALL lbc_lnk( zbathy, 'T', 1. ) 1280 mbathy(:,:) = INT( zbathy(:,:) )1283 mbathy(:,:) = INT( zbathy(:,:) ) 1281 1284 ENDIF 1282 1285 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 1283 misfdep( 1 ,:) = misfdep(jpim1,:) ! local domain is cyclic east-west1286 misfdep( 1 ,:) = misfdep(jpim1,:) ! local domain is cyclic east-west 1284 1287 misfdep(jpi,:) = misfdep( 2 ,:) 1285 ENDIF 1286 1287 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 1288 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west 1289 mbathy(jpi,:) = mbathy( 2 ,:) 1288 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west 1289 mbathy(jpi,:) = mbathy( 2 ,:) 1290 1290 ENDIF 1291 1291 … … 1314 1314 ! find the minimum change option: 1315 1315 ! test bathy 1316 IF (risfdep(ji,jj) .GT.1) THEN1316 IF (risfdep(ji,jj) > 1) THEN 1317 1317 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 1318 1318 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) … … 1320 1320 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1321 1321 1322 IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT.misfdep(ji,jj)) THEN1323 IF (zbathydiff .LE.zrisfdepdiff) THEN1322 IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN 1323 IF (zbathydiff <= zrisfdepdiff) THEN 1324 1324 bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 1325 1325 mbathy(ji,jj)= mbathy(ji,jj) + 1 … … 1338 1338 ! find the minimum change option: 1339 1339 ! test bathy 1340 IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT.1) THEN1340 IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1341 1341 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1)& 1342 1342 & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1343 1343 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) & 1344 1344 & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1345 IF (zbathydiff .LE.zrisfdepdiff) THEN1345 IF (zbathydiff <= zrisfdepdiff) THEN 1346 1346 mbathy(ji,jj) = mbathy(ji,jj) + 1 1347 1347 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) … … 1354 1354 END DO 1355 1355 1356 ! point V mbathy(ji,jj) EQmisfdep(ji,jj+1)1356 ! point V mbathy(ji,jj) == misfdep(ji,jj+1) 1357 1357 DO jj = 1, jpjm1 1358 1358 DO ji = 1, jpim1 1359 IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT.1) THEN1359 IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1360 1360 zbathydiff =ABS(bathy(ji,jj ) - (gdepw_1d(mbathy (ji,jj)+1) & 1361 1361 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat ))) 1362 1362 zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1)) & 1363 1363 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 1364 IF (zbathydiff .LE.zrisfdepdiff) THEN1364 IF (zbathydiff <= zrisfdepdiff) THEN 1365 1365 mbathy(ji,jj) = mbathy(ji,jj) + 1 1366 1366 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) & … … 1376 1376 1377 1377 IF( lk_mpp ) THEN 1378 zbathy(:,:) = FLOAT( misfdep(:,:) )1378 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1379 1379 CALL lbc_lnk( zbathy, 'T', 1. ) 1380 1380 misfdep(:,:) = INT( zbathy(:,:) ) 1381 CALL lbc_lnk( risfdep, 'T', 1. ) 1382 CALL lbc_lnk( bathy, 'T', 1. ) 1383 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1381 1382 CALL lbc_lnk( risfdep,'T', 1. ) 1383 CALL lbc_lnk( bathy, 'T', 1. ) 1384 1385 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1384 1386 CALL lbc_lnk( zbathy, 'T', 1. ) 1385 mbathy(:,:) = INT( zbathy(:,:) )1387 mbathy(:,:) = INT( zbathy(:,:) ) 1386 1388 ENDIF 1387 ! point V misdep(ji,jj) EQmbathy(ji,jj+1)1389 ! point V misdep(ji,jj) == mbathy(ji,jj+1) 1388 1390 DO jj = 1, jpjm1 1389 1391 DO ji = 1, jpim1 1390 IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT.1) THEN1392 IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) > 1) THEN 1391 1393 zbathydiff =ABS( bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) & 1392 1394 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 1393 1395 zrisfdepdiff=ABS(risfdep(ji,jj ) - (gdepw_1d(misfdep(ji,jj ) ) & 1394 1396 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat ))) 1395 IF (zbathydiff .LE.zrisfdepdiff) THEN1397 IF (zbathydiff <= zrisfdepdiff) THEN 1396 1398 mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 1397 1399 bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) … … 1406 1408 1407 1409 IF( lk_mpp ) THEN 1408 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(:,:) ) 1410 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1414 1411 CALL lbc_lnk( zbathy, 'T', 1. ) 1415 mbathy(:,:) = INT( zbathy(:,:) ) 1412 misfdep(:,:) = INT( zbathy(:,:) ) 1413 1414 CALL lbc_lnk( risfdep,'T', 1. ) 1415 CALL lbc_lnk( bathy, 'T', 1. ) 1416 1417 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1418 CALL lbc_lnk( zbathy, 'T', 1. ) 1419 mbathy(:,:) = INT( zbathy(:,:) ) 1416 1420 ENDIF 1417 1421 1418 ! point U mbathy(ji,jj) EQmisfdep(ji,jj+1)1422 ! point U mbathy(ji,jj) == misfdep(ji,jj+1) 1419 1423 DO jj = 1, jpjm1 1420 1424 DO ji = 1, jpim1 1421 IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT.1) THEN1425 IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1422 1426 zbathydiff =ABS( bathy(ji ,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 1423 1427 & + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat ))) 1424 1428 zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) & 1425 1429 & - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 1426 IF (zbathydiff .LE.zrisfdepdiff) THEN1430 IF (zbathydiff <= zrisfdepdiff) THEN 1427 1431 mbathy(ji,jj) = mbathy(ji,jj) + 1 1428 1432 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) … … 1436 1440 1437 1441 IF( lk_mpp ) THEN 1438 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(:,:) ) 1442 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1444 1443 CALL lbc_lnk( zbathy, 'T', 1. ) 1445 mbathy(:,:) = INT( zbathy(:,:) ) 1444 misfdep(:,:) = INT( zbathy(:,:) ) 1445 1446 CALL lbc_lnk( risfdep,'T', 1. ) 1447 CALL lbc_lnk( bathy, 'T', 1. ) 1448 1449 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1450 CALL lbc_lnk( zbathy, 'T', 1. ) 1451 mbathy(:,:) = INT( zbathy(:,:) ) 1446 1452 ENDIF 1447 1453 1448 ! point U misfdep(ji,jj) EQbathy(ji,jj+1)1454 ! point U misfdep(ji,jj) == bathy(ji,jj+1) 1449 1455 DO jj = 1, jpjm1 1450 1456 DO ji = 1, jpim1 1451 IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT.1) THEN1457 IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) > 1) THEN 1452 1458 zbathydiff =ABS( bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) & 1453 1459 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 1454 1460 zrisfdepdiff=ABS(risfdep(ji ,jj) - (gdepw_1d(misfdep(ji ,jj) ) & 1455 1461 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat ))) 1456 IF (zbathydiff .LE.zrisfdepdiff) THEN1462 IF (zbathydiff <= zrisfdepdiff) THEN 1457 1463 mbathy(ji+1,jj) = mbathy (ji+1,jj) + 1 1458 1464 bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) & … … 1468 1474 1469 1475 IF( lk_mpp ) THEN 1470 zbathy(:,:) = FLOAT( misfdep(:,:) )1476 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1471 1477 CALL lbc_lnk( zbathy, 'T', 1. ) 1472 1478 misfdep(:,:) = INT( zbathy(:,:) ) 1473 CALL lbc_lnk( risfdep, 'T', 1. ) 1474 CALL lbc_lnk( bathy, 'T', 1. ) 1475 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1479 1480 CALL lbc_lnk( risfdep,'T', 1. ) 1481 CALL lbc_lnk( bathy, 'T', 1. ) 1482 1483 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1476 1484 CALL lbc_lnk( zbathy, 'T', 1. ) 1477 mbathy(:,:) = INT( zbathy(:,:) )1485 mbathy(:,:) = INT( zbathy(:,:) ) 1478 1486 ENDIF 1479 1487 END DO … … 1485 1493 DO jk = 2, jpk 1486 1494 WHERE (misfdep==0) misfdep=jpk 1487 zmask=0 1488 WHERE (misfdep .LE.jk) zmask=11495 zmask=0._wp 1496 WHERE (misfdep <= jk) zmask=1 1489 1497 DO jj = 2, jpjm1 1490 1498 DO ji = 2, jpim1 1491 IF (misfdep(ji,jj) .EQ.jk) THEN1499 IF (misfdep(ji,jj) == jk) THEN 1492 1500 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1493 IF (ibtest .LE.1) THEN1501 IF (ibtest <= 1) THEN 1494 1502 risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1 1495 IF (misfdep(ji,jj) .GT.mbathy(ji,jj)) misfdep(ji,jj) = jpk1503 IF (misfdep(ji,jj) > mbathy(ji,jj)) misfdep(ji,jj) = jpk 1496 1504 END IF 1497 1505 END IF … … 1500 1508 END DO 1501 1509 WHERE (misfdep==jpk) 1502 misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.1510 misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 1503 1511 END WHERE 1504 1512 IF( lk_mpp ) THEN 1505 zbathy(:,:) = FLOAT( misfdep(:,:) )1513 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1506 1514 CALL lbc_lnk( zbathy, 'T', 1. ) 1507 1515 misfdep(:,:) = INT( zbathy(:,:) ) 1508 CALL lbc_lnk( risfdep, 'T', 1. ) 1509 CALL lbc_lnk( bathy, 'T', 1. ) 1510 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1516 1517 CALL lbc_lnk( risfdep,'T', 1. ) 1518 CALL lbc_lnk( bathy, 'T', 1. ) 1519 1520 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1511 1521 CALL lbc_lnk( zbathy, 'T', 1. ) 1512 mbathy(:,:) = INT( zbathy(:,:) )1522 mbathy(:,:) = INT( zbathy(:,:) ) 1513 1523 ENDIF 1514 1524 1515 1525 ! remove single point "bay" on bathy coast line beneath an ice shelf' 1516 1526 DO jk = jpk,1,-1 1517 zmask=0 1518 WHERE (mbathy .GE.jk ) zmask=11527 zmask=0._wp 1528 WHERE (mbathy >= jk ) zmask=1 1519 1529 DO jj = 2, jpjm1 1520 1530 DO ji = 2, jpim1 1521 IF (mbathy(ji,jj) .EQ. jk .AND. misfdep(ji,jj) .GE.2) THEN1531 IF (mbathy(ji,jj) == jk .AND. misfdep(ji,jj) >= 2) THEN 1522 1532 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1523 IF (ibtest .LE.1) THEN1533 IF (ibtest <= 1) THEN 1524 1534 bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1 1525 IF (misfdep(ji,jj) .GT.mbathy(ji,jj)) mbathy(ji,jj) = 01535 IF (misfdep(ji,jj) > mbathy(ji,jj)) mbathy(ji,jj) = 0 1526 1536 END IF 1527 1537 END IF … … 1530 1540 END DO 1531 1541 WHERE (mbathy==0) 1532 misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.1542 misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 1533 1543 END WHERE 1534 1544 IF( lk_mpp ) THEN 1535 zbathy(:,:) = FLOAT( misfdep(:,:) )1545 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1536 1546 CALL lbc_lnk( zbathy, 'T', 1. ) 1537 1547 misfdep(:,:) = INT( zbathy(:,:) ) 1538 CALL lbc_lnk( risfdep, 'T', 1. ) 1539 CALL lbc_lnk( bathy, 'T', 1. ) 1540 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1548 1549 CALL lbc_lnk( risfdep,'T', 1. ) 1550 CALL lbc_lnk( bathy, 'T', 1. ) 1551 1552 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1541 1553 CALL lbc_lnk( zbathy, 'T', 1. ) 1542 mbathy(:,:) = INT( zbathy(:,:) )1554 mbathy(:,:) = INT( zbathy(:,:) ) 1543 1555 ENDIF 1544 1556 … … 1546 1558 zmisfdep = misfdep 1547 1559 zrisfdep = risfdep 1548 WHERE (zmisfdep .LE. 1) zmisfdep=jpk1560 WHERE (zmisfdep <= 1._wp) zmisfdep=jpk 1549 1561 DO jj = 2, jpjm1 1550 1562 DO ji = 2, jpim1 1551 1563 ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj ) 1552 1564 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)1565 IF( zmisfdep(ji,jj) >= mbathy(ji-1,jj ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj ) - 1) 1566 IF( zmisfdep(ji,jj) >= mbathy(ji+1,jj ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj ) - 1) 1567 IF( zmisfdep(ji,jj) >= mbathy(ji ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji ,jj-1) - 1) 1568 IF( zmisfdep(ji,jj) >= mbathy(ji ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji ,jj+1) - 1) 1557 1569 ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1558 IF( ibtest == jpk .AND. misfdep(ji,jj) .GE.2) THEN1570 IF( ibtest == jpk .AND. misfdep(ji,jj) >= 2) THEN 1559 1571 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp 1560 1572 END IF 1561 IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) .GE.2) THEN1573 IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) >= 2) THEN 1562 1574 misfdep(ji,jj) = ibtest 1563 1575 risfdep(ji,jj) = gdepw_1d(ibtest) … … 1567 1579 1568 1580 IF( lk_mpp ) THEN 1569 zbathy(:,:) = FLOAT( misfdep(:,:) )1570 CALL lbc_lnk( zbathy, 'T', 1. )1581 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1582 CALL lbc_lnk( zbathy, 'T', 1. ) 1571 1583 misfdep(:,:) = INT( zbathy(:,:) ) 1584 1572 1585 CALL lbc_lnk( risfdep, 'T', 1. ) 1573 CALL lbc_lnk( bathy, 'T', 1. ) 1586 CALL lbc_lnk( bathy, 'T', 1. ) 1587 1574 1588 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1575 CALL lbc_lnk( zbathy, 'T', 1. )1589 CALL lbc_lnk( zbathy, 'T', 1. ) 1576 1590 mbathy(:,:) = INT( zbathy(:,:) ) 1577 1591 ENDIF … … 1583 1597 ibtestim1 = zmbathy(ji-1,jj ) ; ibtestip1 = zmbathy(ji+1,jj ) 1584 1598 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 = 01599 IF( zmbathy(ji,jj) < misfdep(ji-1,jj ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj ) + 1) 1600 IF( zmbathy(ji,jj) < misfdep(ji+1,jj ) ) ibtestip1 = 0 1601 IF( zmbathy(ji,jj) < misfdep(ji ,jj-1) ) ibtestjm1 = 0 1602 IF( zmbathy(ji,jj) < misfdep(ji ,jj+1) ) ibtestjp1 = 0 1589 1603 ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1590 IF( ibtest == 0 .AND. misfdep(ji,jj) .GE.2) THEN1604 IF( ibtest == 0 .AND. misfdep(ji,jj) >= 2) THEN 1591 1605 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 1592 1606 END IF 1593 IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) .GE.2) THEN1607 IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) >= 2) THEN 1594 1608 mbathy(ji,jj) = ibtest 1595 1609 bathy(ji,jj) = gdepw_1d(ibtest+1) … … 1598 1612 END DO 1599 1613 IF( lk_mpp ) THEN 1600 zbathy(:,:) = FLOAT( misfdep(:,:) )1601 CALL lbc_lnk( zbathy, 'T', 1. )1614 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1615 CALL lbc_lnk( zbathy, 'T', 1. ) 1602 1616 misfdep(:,:) = INT( zbathy(:,:) ) 1617 1603 1618 CALL lbc_lnk( risfdep, 'T', 1. ) 1604 CALL lbc_lnk( bathy, 'T', 1. ) 1619 CALL lbc_lnk( bathy, 'T', 1. ) 1620 1605 1621 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1606 CALL lbc_lnk( zbathy, 'T', 1. )1622 CALL lbc_lnk( zbathy, 'T', 1. ) 1607 1623 mbathy(:,:) = INT( zbathy(:,:) ) 1608 1624 ENDIF … … 1610 1626 DO jj = 1, jpjm1 1611 1627 DO ji = 1, jpim1 1612 IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE.1) THEN1628 IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 1613 1629 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1614 1630 END IF … … 1616 1632 END DO 1617 1633 IF( lk_mpp ) THEN 1618 zbathy(:,:) = FLOAT( misfdep(:,:) )1619 CALL lbc_lnk( zbathy, 'T', 1. )1634 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1635 CALL lbc_lnk( zbathy, 'T', 1. ) 1620 1636 misfdep(:,:) = INT( zbathy(:,:) ) 1637 1621 1638 CALL lbc_lnk( risfdep, 'T', 1. ) 1622 CALL lbc_lnk( bathy, 'T', 1. ) 1639 CALL lbc_lnk( bathy, 'T', 1. ) 1640 1623 1641 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1624 CALL lbc_lnk( zbathy, 'T', 1. )1642 CALL lbc_lnk( zbathy, 'T', 1. ) 1625 1643 mbathy(:,:) = INT( zbathy(:,:) ) 1626 1644 ENDIF … … 1628 1646 DO jj = 1, jpjm1 1629 1647 DO ji = 1, jpim1 1630 IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE.1) THEN1648 IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 1631 1649 mbathy(ji+1,jj) = mbathy(ji+1,jj) - 1; bathy(ji+1,jj) = gdepw_1d(mbathy(ji+1,jj)+1) ; 1632 1650 END IF … … 1634 1652 END DO 1635 1653 IF( lk_mpp ) THEN 1636 zbathy(:,:) = FLOAT( misfdep(:,:) )1654 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1637 1655 CALL lbc_lnk( zbathy, 'T', 1. ) 1638 1656 misfdep(:,:) = INT( zbathy(:,:) ) 1639 CALL lbc_lnk( risfdep, 'T', 1. ) 1640 CALL lbc_lnk( bathy, 'T', 1. ) 1657 1658 CALL lbc_lnk( risfdep,'T', 1. ) 1659 CALL lbc_lnk( bathy, 'T', 1. ) 1660 1641 1661 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1642 1662 CALL lbc_lnk( zbathy, 'T', 1. ) … … 1646 1666 DO jj = 1, jpjm1 1647 1667 DO ji = 1, jpi 1648 IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE.1) THEN1668 IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 1649 1669 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1650 1670 END IF … … 1652 1672 END DO 1653 1673 IF( lk_mpp ) THEN 1654 zbathy(:,:) = FLOAT( misfdep(:,:) )1674 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1655 1675 CALL lbc_lnk( zbathy, 'T', 1. ) 1656 1676 misfdep(:,:) = INT( zbathy(:,:) ) 1657 CALL lbc_lnk( risfdep, 'T', 1. ) 1658 CALL lbc_lnk( bathy, 'T', 1. ) 1677 1678 CALL lbc_lnk( risfdep,'T', 1. ) 1679 CALL lbc_lnk( bathy, 'T', 1. ) 1680 1659 1681 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1660 1682 CALL lbc_lnk( zbathy, 'T', 1. ) … … 1664 1686 DO jj = 1, jpjm1 1665 1687 DO ji = 1, jpi 1666 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) 1688 IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 1689 mbathy(ji,jj+1) = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ; 1668 1690 END IF 1669 1691 END DO 1670 1692 END DO 1671 1693 IF( lk_mpp ) THEN 1672 zbathy(:,:) = FLOAT( misfdep(:,:) )1694 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1673 1695 CALL lbc_lnk( zbathy, 'T', 1. ) 1674 1696 misfdep(:,:) = INT( zbathy(:,:) ) 1675 CALL lbc_lnk( risfdep, 'T', 1. ) 1676 CALL lbc_lnk( bathy, 'T', 1. ) 1697 1698 CALL lbc_lnk( risfdep,'T', 1. ) 1699 CALL lbc_lnk( bathy, 'T', 1. ) 1700 1677 1701 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1678 1702 CALL lbc_lnk( zbathy, 'T', 1. ) … … 1694 1718 ! end check compatibility ice shelf/bathy 1695 1719 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1696 WHERE ( misfdep(:,:) <= 5)1720 WHERE (risfdep(:,:) <= 10._wp) 1697 1721 misfdep = 1; risfdep = 0.0_wp; 1698 1722 END WHERE … … 1814 1838 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf') 1815 1839 1816 END SUBROUTINE 1840 END SUBROUTINE zgr_isf 1817 1841 1818 1842 SUBROUTINE zgr_sco … … 2143 2167 fsde3w(:,:,:) = gdep3w_0(:,:,:) 2144 2168 ! 2145 where (e3t_0 (:,:,:) .eq.0.0) e3t_0(:,:,:) = 1.02146 where (e3u_0 (:,:,:) .eq.0.0) e3u_0(:,:,:) = 1.02147 where (e3v_0 (:,:,:) .eq.0.0) e3v_0(:,:,:) = 1.02148 where (e3f_0 (:,:,:) .eq.0.0) e3f_0(:,:,:) = 1.02149 where (e3w_0 (:,:,:) .eq.0.0) e3w_0(:,:,:) = 1.02150 where (e3uw_0 (:,:,:) .eq.0.0) e3uw_0(:,:,:) = 1.02151 where (e3vw_0 (:,:,:) .eq.0.0) e3vw_0(:,:,:) = 1.02169 where (e3t_0 (:,:,:) == 0.0) e3t_0(:,:,:) = 1.0_wp 2170 where (e3u_0 (:,:,:) == 0.0) e3u_0(:,:,:) = 1.0_wp 2171 where (e3v_0 (:,:,:) == 0.0) e3v_0(:,:,:) = 1.0_wp 2172 where (e3f_0 (:,:,:) == 0.0) e3f_0(:,:,:) = 1.0_wp 2173 where (e3w_0 (:,:,:) == 0.0) e3w_0(:,:,:) = 1.0_wp 2174 where (e3uw_0 (:,:,:) == 0.0) e3uw_0(:,:,:) = 1.0_wp 2175 where (e3vw_0 (:,:,:) == 0.0) e3vw_0(:,:,:) = 1.0_wp 2152 2176 2153 2177 #if defined key_agrif -
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r5921 r5944 235 235 & * (1._wp - tmask(ji,jj,jk)) 236 236 END DO 237 IF (ikt .GE.2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) &237 IF (ikt >= 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) & 238 238 & * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 239 239 END DO … … 521 521 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 522 522 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztstop 523 REAL(wp), POINTER, DIMENSION(:,:) :: zrhdtop_oce , zpshpi, zpshpj523 REAL(wp), POINTER, DIMENSION(:,:) :: zrhdtop_oce 524 524 !!---------------------------------------------------------------------- 525 525 ! 526 526 CALL wrk_alloc( jpi,jpj, 2, ztstop) 527 527 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj) 528 CALL wrk_alloc( jpi,jpj, zrhdtop_oce , zpshpi, zpshpj)528 CALL wrk_alloc( jpi,jpj, zrhdtop_oce ) 529 529 ! 530 530 ! Local constant initialization … … 605 605 CALL wrk_dealloc( jpi,jpj,2 , ztstop) 606 606 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj) 607 CALL wrk_dealloc( jpi,jpj , zrhdtop_oce , zpshpi, zpshpj)607 CALL wrk_dealloc( jpi,jpj , zrhdtop_oce ) 608 608 ! 609 609 END SUBROUTINE hpg_isf -
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r5624 r5944 100 100 ! 101 101 INTEGER :: ji, jj, jk ! dummy loop indices 102 INTEGER :: ik u, ikv! local integers102 INTEGER :: ikt ! local integers 103 103 #if ! defined key_dynspg_flt 104 104 REAL(wp) :: z2dt ! temporary scalar … … 268 268 DO jj = 1,jpj 269 269 DO ji = 1,jpi 270 jk= mikt(ji,jj)271 fse3t_b(ji,jj, jk) = fse3t_b(ji,jj,jk) - atfp * rdt * r1_rau0 &270 ikt = mikt(ji,jj) 271 fse3t_b(ji,jj,ikt) = fse3t_b(ji,jj,ikt) - atfp * rdt * r1_rau0 & 272 272 & * ( (emp_b(ji,jj) - emp(ji,jj) ) & 273 273 & - (rnf_b(ji,jj) - rnf(ji,jj) ) & 274 & + (fwfisf_b(ji,jj) - fwfisf(ji,jj) ) ) * tmask(ji,jj, jk)274 & + (fwfisf_b(ji,jj) - fwfisf(ji,jj) ) ) * tmask(ji,jj,ikt) 275 275 END DO 276 276 END DO -
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r5910 r5944 934 934 ! 935 935 IF ( (.NOT.Agrif_Root()).AND.(ln_bt_fw) ) THEN 936 IF ( Agrif_NbStepint() .EQ.0 ) THEN936 IF ( Agrif_NbStepint() == 0 ) THEN 937 937 ub2_i_b(:,:) = 0.e0 938 938 vb2_i_b(:,:) = 0.e0 -
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90
r5189 r5944 86 86 DO jj = 2, jpj ! vertical fluxes 87 87 DO ji = fs_2, jpi 88 zww(ji,jj) = 0.25_wp * e1 e2t(ji,jj) * wn(ji,jj,jk)88 zww(ji,jj) = 0.25_wp * e12t(ji,jj) * wn(ji,jj,jk) 89 89 END DO 90 90 END DO … … 193 193 DO jj = 2, jpj 194 194 DO ji = fs_2, jpi ! vector opt. 195 zww(ji,jj,jk) = 0.25_wp * e1 e2t(ji,jj) * wn(ji,jj,jk)195 zww(ji,jj,jk) = 0.25_wp * e12t(ji,jj) * wn(ji,jj,jk) 196 196 END DO 197 197 END DO -
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
r5200 r5944 586 586 ! 587 587 jk = nmln(ji+ip,jj) + 1 588 IF( jk .GT.mbkt(ji+ip,jj) ) THEN !ML reaches bottom588 IF( jk > mbkt(ji+ip,jj) ) THEN !ML reaches bottom 589 589 zti_mlb(ji+ip,jj ,1-ip,kp) = 0.0_wp 590 590 ELSE … … 596 596 ! 597 597 jk = nmln(ji,jj+jp) + 1 598 IF( jk .GT.mbkt(ji,jj+jp) ) THEN !ML reaches bottom598 IF( jk > mbkt(ji,jj+jp) ) THEN !ML reaches bottom 599 599 ztj_mlb(ji ,jj+jp,1-jp,kp) = 0.0_wp 600 600 ELSE -
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90
r5624 r5944 91 91 IF( kn_fwb == 3 .AND. ln_isfcav ) CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 with ln_isfcav = .TRUE. not working, we stop ' ) 92 92 ! 93 area = glob_sum( e1 e2t(:,:) * tmask(:,:,1)) ! interior global domain surface93 area = glob_sum( e12t(:,:) * tmask(:,:,1)) ! interior global domain surface 94 94 ! isf cavities are excluded because it can feedback to the melting with generation of inhibition of plumes 95 95 ! and in case of no melt, it can generate HSSW. … … 108 108 ! 109 109 IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN 110 z_fwf = glob_sum( e1 e2t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) - snwice_fmass(:,:) ) ) / area ! sum over the global domain110 z_fwf = glob_sum( e12t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) - snwice_fmass(:,:) ) ) / area ! sum over the global domain 111 111 zcoef = z_fwf * rcp 112 112 emp(:,:) = emp(:,:) - z_fwf * tmask(:,:,1) … … 133 133 a_fwb_b = a_fwb ! mean sea level taking into account the ice+snow 134 134 ! sum over the global domain 135 a_fwb = glob_sum( e1 e2t(:,:) * ( sshn(:,:) + snwice_mass(:,:) * r1_rau0 ) )135 a_fwb = glob_sum( e12t(:,:) * ( sshn(:,:) + snwice_mass(:,:) * r1_rau0 ) ) 136 136 a_fwb = a_fwb * 1.e+3 / ( area * rday * 365. ) ! convert in Kg/m3/s = mm/s 137 137 !!gm ! !!bug 365d year … … 159 159 ztmsk_neg(:,:) = tmask_i(:,:) - ztmsk_pos(:,:) 160 160 ! 161 zsurf_neg = glob_sum( e1 e2t(:,:)*ztmsk_neg(:,:) ) ! Area filled by <0 and >0 erp162 zsurf_pos = glob_sum( e1 e2t(:,:)*ztmsk_pos(:,:) )161 zsurf_neg = glob_sum( e12t(:,:)*ztmsk_neg(:,:) ) ! Area filled by <0 and >0 erp 162 zsurf_pos = glob_sum( e12t(:,:)*ztmsk_pos(:,:) ) 163 163 ! ! fwf global mean (excluding ocean to ice/snow exchanges) 164 z_fwf = glob_sum( e1 e2t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) - snwice_fmass(:,:) ) ) / area164 z_fwf = glob_sum( e12t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) - snwice_fmass(:,:) ) ) / area 165 165 ! 166 166 IF( z_fwf < 0._wp ) THEN ! spread out over >0 erp area to increase evaporation … … 172 172 ENDIF 173 173 ! 174 zsum_fwf = glob_sum( e1 e2t(:,:) * z_fwf ) ! fwf global mean over <0 or >0 erp area174 zsum_fwf = glob_sum( e12t(:,:) * z_fwf ) ! fwf global mean over <0 or >0 erp area 175 175 !!gm : zsum_fwf = z_fwf * area ??? it is right? I think so.... 176 176 z_fwf_nsrf = zsum_fwf / ( zsurf_tospread + rsmall ) 177 177 ! ! weight to respect erp field 2D structure 178 zsum_erp = glob_sum( ztmsk_tospread(:,:) * erp(:,:) * e1 e2t(:,:) )178 zsum_erp = glob_sum( ztmsk_tospread(:,:) * erp(:,:) * e12t(:,:) ) 179 179 z_wgt(:,:) = ztmsk_tospread(:,:) * erp(:,:) / ( zsum_erp + rsmall ) 180 180 ! ! final correction term to apply … … 191 191 IF( z_fwf < 0._wp ) THEN 192 192 WRITE(numout,*)' z_fwf < 0' 193 WRITE(numout,*)' SUM(erp+) = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1 e2t(:,:) )*1.e-9,' Sv'193 WRITE(numout,*)' SUM(erp+) = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e12t(:,:) )*1.e-9,' Sv' 194 194 ELSE 195 195 WRITE(numout,*)' z_fwf >= 0' 196 WRITE(numout,*)' SUM(erp-) = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1 e2t(:,:) )*1.e-9,' Sv'196 WRITE(numout,*)' SUM(erp-) = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e12t(:,:) )*1.e-9,' Sv' 197 197 ENDIF 198 WRITE(numout,*)' SUM(empG) = ', SUM( z_fwf*e1 e2t(:,:) )*1.e-9,' Sv'198 WRITE(numout,*)' SUM(empG) = ', SUM( z_fwf*e12t(:,:) )*1.e-9,' Sv' 199 199 WRITE(numout,*)' z_fwf = ', z_fwf ,' Kg/m2/s' 200 200 WRITE(numout,*)' z_fwf_nsrf = ', z_fwf_nsrf ,' Kg/m2/s' -
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
r5905 r5944 56 56 !: (first wet level and last level include in the tbl) 57 57 #else 58 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: misfkt, misfkb !:Level of ice shelf base58 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: misfkt, misfkb !:Level of ice shelf base 59 59 #endif 60 60 61 61 62 REAL(wp), PUBLIC, SAVE :: rcpi = 2000.0_wp ! phycst ?63 REAL(wp), PUBLIC, SAVE :: kappa= 1.54e-6_wp ! phycst ?64 REAL(wp), PUBLIC, SAVE :: rhoisf = 920.0_wp ! phycst ?65 REAL(wp), PUBLIC, SAVE :: tsurf = -20.0_wp ! phycst ?66 REAL(wp), PUBLIC, SAVE :: lfusisf= 0.334e6_wp ! phycst ?62 REAL(wp), PUBLIC, SAVE :: rcpi = 2000.0_wp ! phycst ? 63 REAL(wp), PUBLIC, SAVE :: rkappa = 1.54e-6_wp ! phycst ? 64 REAL(wp), PUBLIC, SAVE :: rhoisf = 920.0_wp ! phycst ? 65 REAL(wp), PUBLIC, SAVE :: tsurf = -20.0_wp ! phycst ? 66 REAL(wp), PUBLIC, SAVE :: rlfusisf = 0.334e6_wp ! phycst ? 67 67 68 68 !: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3) 69 CHARACTER(len=100), PUBLIC :: cn_dirisf = './' !: Root directory for location of ssr files 70 TYPE(FLD_N) , PUBLIC :: sn_fwfisf !: information about the isf file to be read 71 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_fwfisf 72 TYPE(FLD_N) , PUBLIC :: sn_rnfisf !: information about the runoff file to be read 73 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnfisf 74 TYPE(FLD_N) , PUBLIC :: sn_depmax_isf, sn_depmin_isf, sn_Leff_isf !: information about the runoff file to be read 69 CHARACTER(len=100), PUBLIC :: cn_dirisf = './' !: Root directory for location of ssr files 70 TYPE(FLD_N) , PUBLIC :: sn_fwfisf !: information about the isf file to be read 71 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_fwfisf 72 TYPE(FLD_N) , PUBLIC :: sn_rnfisf !: information about the runoff file to be read 73 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnfisf 74 TYPE(FLD_N) , PUBLIC :: sn_depmax_isf !: information about the runoff file to be read ?? 75 TYPE(FLD_N) , PUBLIC :: sn_depmin_isf !: information about the runoff file to be read ?? 76 TYPE(FLD_N) , PUBLIC :: sn_Leff_isf !: information about the runoff file to be read ?? 75 77 76 78 !! * Substitutions … … 85 87 86 88 SUBROUTINE sbc_isf(kt) 87 INTEGER, INTENT(in) :: kt ! ocean time step88 INTEGER :: ji, jj, jk ! loop index89 INTEGER :: ikt, ikb ! top and bottom level of the isf boundary layer90 INTEGER :: inum, ierror91 INTEGER :: ios ! Local integer output status for namelist read92 REAL(wp) :: zhk93 REAL(wp), DIMENSION (:,:), POINTER :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep)94 CHARACTER(len=256) :: cvarzisf, cvarhisf ! name for isf file95 CHARACTER(LEN=32 ) :: cvarLeff ! variable name for efficient Length scale96 !97 89 !!--------------------------------------------------------------------- 98 NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, nn_isf , & 99 & sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf 100 ! 101 ! allocation 102 CALL wrk_alloc( jpi,jpj, zt_frz, zdep ) 90 !! *** ROUTINE sbc_isf *** 91 !! 92 !! ** Purpose : Compute Salt and Heat fluxes related to ice_shelf 93 !! melting and freezing 94 !! 95 !! ** Method : 4 parameterizations are available according to nn_isf 96 !! nn_isf = 1 : Realistic ice_shelf formulation 97 !! 2 : Beckmann & Goose parameterization 98 !! 3 : Specified runoff in deptht (Mathiot & al. ) 99 !! 4 : specified fwf and heat flux forcing beneath the ice shelf 100 !!---------------------------------------------------------------------- 101 INTEGER, INTENT( in ) :: kt ! ocean time step 102 ! 103 INTEGER :: ji, jj ! loop index 104 REAL(wp), DIMENSION (:,:), POINTER :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep) 105 !!--------------------------------------------------------------------- 103 106 ! 104 107 ! ! ====================== ! 105 108 IF( kt == nit000 ) THEN ! First call kt=nit000 ! 106 109 ! ! ====================== ! 107 REWIND( numnam_ref ) ! Namelist namsbc_rnf in reference namelist : Runoffs 108 READ ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901) 109 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp ) 110 111 REWIND( numnam_cfg ) ! Namelist namsbc_rnf in configuration namelist : Runoffs 112 READ ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 ) 113 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp ) 114 IF(lwm) WRITE ( numond, namsbc_isf ) 115 116 117 IF ( lwp ) WRITE(numout,*) 118 IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf' 119 IF ( lwp ) WRITE(numout,*) '~~~~~~~~~' 120 IF ( lwp ) WRITE(numout,*) 'sbcisf :' 121 IF ( lwp ) WRITE(numout,*) '~~~~~~~~' 122 IF ( lwp ) WRITE(numout,*) ' nn_isf = ', nn_isf 123 IF ( lwp ) WRITE(numout,*) ' nn_isfblk = ', nn_isfblk 124 IF ( lwp ) WRITE(numout,*) ' rn_hisf_tbl = ', rn_hisf_tbl 125 IF ( lwp ) WRITE(numout,*) ' nn_gammablk = ', nn_gammablk 126 IF ( lwp ) WRITE(numout,*) ' rn_gammat0 = ', rn_gammat0 127 IF ( lwp ) WRITE(numout,*) ' rn_gammas0 = ', rn_gammas0 128 IF ( lwp ) WRITE(numout,*) ' rn_tfri2 = ', rn_tfri2 110 CALL sbc_isf_init 111 ! ! ---------------------------------------- ! 112 ELSE ! Swap of forcing fields ! 113 ! ! ---------------------------------------- ! 114 fwfisf_b (:,: ) = fwfisf (:,: ) ! Swap the ocean forcing fields except at nit000 115 risf_tsc_b(:,:,:) = risf_tsc(:,:,:) ! where before fields are set at the end of the routine 129 116 ! 130 ! Allocate public variable131 IF ( sbc_isf_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' )132 !133 ! initialisation134 qisf(:,:) = 0._wp ; fwfisf (:,:) = 0._wp135 risf_tsc(:,:,:) = 0._wp ; fwfisf_b(:,:) = 0._wp136 !137 ! define isf tbl tickness, top and bottom indice138 IF (nn_isf == 1) THEN139 rhisf_tbl(:,:) = rn_hisf_tbl140 misfkt(:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv141 ELSE IF ((nn_isf == 3) .OR. (nn_isf == 2)) THEN142 ALLOCATE( sf_rnfisf(1), STAT=ierror )143 ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) )144 CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )145 146 !: read effective lenght (BG03)147 IF (nn_isf == 2) THEN148 cvarLeff = 'soLeff'149 CALL iom_open( sn_Leff_isf%clname, inum )150 CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1)151 CALL iom_close(inum)152 !153 risfLeff = risfLeff*1000.0_wp !: convertion in m154 END IF155 156 ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth)157 CALL iom_open( sn_depmax_isf%clname, inum )158 cvarhisf = TRIM(sn_depmax_isf%clvar)159 CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base160 CALL iom_close(inum)161 !162 CALL iom_open( sn_depmin_isf%clname, inum )163 cvarzisf = TRIM(sn_depmin_isf%clvar)164 CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base165 CALL iom_close(inum)166 !167 rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:) !: tickness isf boundary layer168 169 !! compute first level of the top boundary layer170 DO ji = 1, jpi171 DO jj = 1, jpj172 jk = 2173 DO WHILE ( jk .LE. mbkt(ji,jj) .AND. fsdepw(ji,jj,jk) < rzisf_tbl(ji,jj) ) ; jk = jk + 1 ; END DO174 misfkt(ji,jj) = jk-1175 END DO176 END DO177 178 ELSE IF ( nn_isf == 4 ) THEN179 ! as in nn_isf == 1180 rhisf_tbl(:,:) = rn_hisf_tbl181 misfkt(:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv182 183 ! load variable used in fldread (use for temporal interpolation of isf fwf forcing)184 ALLOCATE( sf_fwfisf(1), STAT=ierror )185 ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) )186 CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )187 END IF188 189 rhisf_tbl_0(:,:) = rhisf_tbl(:,:)190 191 ! compute bottom level of isf tbl and thickness of tbl below the ice shelf192 DO jj = 1,jpj193 DO ji = 1,jpi194 ikt = misfkt(ji,jj)195 ikb = misfkt(ji,jj)196 ! thickness of boundary layer at least the top level thickness197 rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t_n(ji,jj,ikt))198 199 ! determine the deepest level influenced by the boundary layer200 DO jk = ikt+1, mbkt(ji,jj)201 IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk202 END DO203 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.204 misfkb(ji,jj) = ikb ! last wet level of the tbl205 r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)206 207 zhk = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1208 ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb) ! proportion of bottom cell influenced by boundary layer209 END DO210 END DO211 212 117 END IF 213 118 214 ! ! ---------------------------------------- !215 IF( kt /= nit000 ) THEN ! Swap of forcing fields !216 ! ! ---------------------------------------- !217 fwfisf_b (:,: ) = fwfisf (:,: ) ! Swap the ocean forcing fields except at nit000218 risf_tsc_b(:,:,:) = risf_tsc(:,:,:) ! where before fields are set at the end of the routine219 !220 ENDIF221 222 119 IF( MOD( kt-1, nn_fsbc) == 0 ) THEN 223 224 225 ! compute salf and heat flux 226 IF (nn_isf == 1) THEN 227 ! realistic ice shelf formulation 120 ! allocation 121 CALL wrk_alloc( jpi,jpj, zt_frz, zdep ) 122 123 ! compute salt and heat flux 124 SELECT CASE ( nn_isf ) 125 CASE ( 1 ) ! realistic ice shelf formulation 228 126 ! compute T/S/U/V for the top boundary layer 229 127 CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T') … … 239 137 CALL sbc_isf_cav (kt) 240 138 241 ELSE IF (nn_isf == 2) THEN 242 ! Beckmann and Goosse parametrisation 139 CASE ( 2 ) ! Beckmann and Goosse parametrisation 243 140 stbl(:,:) = soce 244 141 CALL sbc_isf_bg03(kt) 245 142 246 ELSE IF (nn_isf == 3) THEN 247 ! specified runoff in depth (Mathiot et al., XXXX in preparation) 143 CASE ( 3 ) ! specified runoff in depth (Mathiot et al., XXXX in preparation) 248 144 CALL fld_read ( kt, nn_fsbc, sf_rnfisf ) 249 145 fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1) ! fwf flux from the isf (fwfisf <0 mean melting) 250 qisf(:,:) = fwfisf(:,:) * lfusisf! heat flux146 qisf(:,:) = fwfisf(:,:) * rlfusisf ! heat flux 251 147 stbl(:,:) = soce 252 148 253 ELSE IF (nn_isf == 4) THEN 254 ! specified fwf and heat flux forcing beneath the ice shelf 149 CASE ( 4 ) ! specified fwf and heat flux forcing beneath the ice shelf 255 150 CALL fld_read ( kt, nn_fsbc, sf_fwfisf ) 256 151 fwfisf(:,:) = - sf_fwfisf(1)%fnow(:,:,1) ! fwf flux from the isf (fwfisf <0 mean melting) 257 qisf(:,:) = fwfisf(:,:) * lfusisf! heat flux152 qisf(:,:) = fwfisf(:,:) * rlfusisf ! heat flux 258 153 stbl(:,:) = soce 259 END IF 154 155 END SELECT 260 156 261 157 ! compute tsc due to isf … … 276 172 CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.) 277 173 CALL lbc_lnk(risf_tsc(:,:,jp_sal),'T',1.) 278 CALL lbc_lnk(fwfisf(:,:) ,'T',1.)279 CALL lbc_lnk(qisf(:,:) ,'T',1.)174 CALL lbc_lnk(fwfisf(:,:) ,'T',1.) 175 CALL lbc_lnk(qisf(:,:) ,'T',1.) 280 176 281 177 IF( kt == nit000 ) THEN ! set the forcing field at nit000 - 1 ! … … 290 186 risf_tsc_b(:,:,:)= risf_tsc(:,:,:) 291 187 END IF 292 END IF188 END IF 293 189 ! 294 190 ! output 295 IF( iom_use('qisf' ) ) CALL iom_put('qisf' , qisf) 296 IF( iom_use('fwfisf') ) CALL iom_put('fwfisf', fwfisf) 191 ! JMM : iom_use not necessary here, qisf and fwfisf are always computed. 192 ! If not required in iodef.xml, iom_put does not do anything 193 ! IF( iom_use('qisf' ) ) CALL iom_put('qisf' , qisf) 194 ! IF( iom_use('fwfisf') ) CALL iom_put('fwfisf', fwfisf) 195 CALL iom_put('qisf' , qisf) 196 CALL iom_put('fwfisf', fwfisf) 197 198 ! deallocation 199 CALL wrk_dealloc( jpi,jpj, zt_frz, zdep ) 297 200 END IF 298 299 ! deallocation300 CALL wrk_dealloc( jpi,jpj, zt_frz, zdep )301 201 302 202 END SUBROUTINE sbc_isf … … 315 215 & STAT= sbc_isf_alloc ) 316 216 ! 317 IF( lk_mpp 217 IF( lk_mpp ) CALL mpp_sum ( sbc_isf_alloc ) 318 218 IF( sbc_isf_alloc /= 0 ) CALL ctl_warn('sbc_isf_alloc: failed to allocate arrays.') 319 219 ! 320 END IF220 END IF 321 221 END FUNCTION 322 222 223 SUBROUTINE sbc_isf_init 224 !!--------------------------------------------------------------------- 225 !! *** ROUTINE sbc_isf_init *** 226 !! 227 !! ** Purpose : Initialisation of variables for iceshelf fluxes formulation 228 !! 229 !! ** Method : 4 parameterizations are available according to nn_isf 230 !! nn_isf = 1 : Realistic ice_shelf formulation 231 !! 2 : Beckmann & Goose parameterization 232 !! 3 : Specified runoff in deptht (Mathiot & al. ) 233 !! 4 : specified fwf and heat flux forcing beneath the ice shelf 234 !!---------------------------------------------------------------------- 235 INTEGER :: ji, jj, jk ! loop index 236 INTEGER :: ik ! current level index 237 INTEGER :: ikt, ikb ! top and bottom level of the isf boundary layer 238 INTEGER :: inum, ierror 239 INTEGER :: ios ! Local integer output status for namelist read 240 REAL(wp) :: zhk 241 CHARACTER(len=256) :: cvarzisf, cvarhisf ! name for isf file 242 CHARACTER(LEN=32 ) :: cvarLeff ! variable name for efficient Length scale 243 !!---------------------------------------------------------------------- 244 NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, nn_isf, & 245 & sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf 246 !!---------------------------------------------------------------------- 247 248 REWIND( numnam_ref ) ! Namelist namsbc_rnf in reference namelist : Runoffs 249 READ ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901) 250 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp ) 251 252 REWIND( numnam_cfg ) ! Namelist namsbc_rnf in configuration namelist : Runoffs 253 READ ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 ) 254 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp ) 255 IF(lwm) WRITE ( numond, namsbc_isf ) 256 257 IF ( lwp ) WRITE(numout,*) 258 IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf' 259 IF ( lwp ) WRITE(numout,*) '~~~~~~~~~' 260 IF ( lwp ) WRITE(numout,*) 'sbcisf :' 261 IF ( lwp ) WRITE(numout,*) '~~~~~~~~' 262 IF ( lwp ) WRITE(numout,*) ' nn_isf = ', nn_isf 263 IF ( lwp ) WRITE(numout,*) ' nn_isfblk = ', nn_isfblk 264 IF ( lwp ) WRITE(numout,*) ' rn_hisf_tbl = ', rn_hisf_tbl 265 IF ( lwp ) WRITE(numout,*) ' nn_gammablk = ', nn_gammablk 266 IF ( lwp ) WRITE(numout,*) ' rn_gammat0 = ', rn_gammat0 267 IF ( lwp ) WRITE(numout,*) ' rn_gammas0 = ', rn_gammas0 268 IF ( lwp ) WRITE(numout,*) ' rn_tfri2 = ', rn_tfri2 269 ! 270 ! Allocate public variable 271 IF ( sbc_isf_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' ) 272 ! 273 ! initialisation 274 qisf(:,:) = 0._wp ; fwfisf (:,:) = 0._wp 275 risf_tsc(:,:,:) = 0._wp ; fwfisf_b(:,:) = 0._wp 276 ! 277 ! define isf tbl tickness, top and bottom indice 278 SELECT CASE ( nn_isf ) 279 CASE ( 1 ) 280 rhisf_tbl(:,:) = rn_hisf_tbl 281 misfkt(:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv 282 283 CASE ( 2 , 3 ) 284 ALLOCATE( sf_rnfisf(1), STAT=ierror ) 285 ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) ) 286 CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 287 288 ! read effective lenght (BG03) 289 IF (nn_isf == 2) THEN 290 CALL iom_open( sn_Leff_isf%clname, inum ) 291 cvarLeff = TRIM(sn_Leff_isf%clvar) 292 CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1) 293 CALL iom_close(inum) 294 ! 295 risfLeff = risfLeff*1000.0_wp !: convertion in m 296 END IF 297 298 ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth) 299 CALL iom_open( sn_depmax_isf%clname, inum ) 300 cvarhisf = TRIM(sn_depmax_isf%clvar) 301 CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base 302 CALL iom_close(inum) 303 ! 304 CALL iom_open( sn_depmin_isf%clname, inum ) 305 cvarzisf = TRIM(sn_depmin_isf%clvar) 306 CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base 307 CALL iom_close(inum) 308 ! 309 rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:) !: tickness isf boundary layer 310 311 !! compute first level of the top boundary layer 312 DO ji = 1, jpi 313 DO jj = 1, jpj 314 ik = 2 315 DO WHILE ( ik <= mbkt(ji,jj) .AND. fsdepw(ji,jj,ik) < rzisf_tbl(ji,jj) ) ; ik = ik + 1 ; END DO 316 misfkt(ji,jj) = ik-1 317 END DO 318 END DO 319 320 CASE ( 4 ) 321 ! as in nn_isf == 1 322 rhisf_tbl(:,:) = rn_hisf_tbl 323 misfkt(:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv 324 325 ! load variable used in fldread (use for temporal interpolation of isf fwf forcing) 326 ALLOCATE( sf_fwfisf(1), STAT=ierror ) 327 ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) ) 328 CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 329 330 END SELECT 331 332 rhisf_tbl_0(:,:) = rhisf_tbl(:,:) 333 334 ! compute bottom level of isf tbl and thickness of tbl below the ice shelf 335 DO jj = 1,jpj 336 DO ji = 1,jpi 337 ikt = misfkt(ji,jj) 338 ikb = misfkt(ji,jj) 339 ! thickness of boundary layer at least the top level thickness 340 rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t_n(ji,jj,ikt)) 341 342 ! determine the deepest level influenced by the boundary layer 343 DO jk = ikt+1, mbkt(ji,jj) 344 IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 345 END DO 346 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 347 misfkb(ji,jj) = ikb ! last wet level of the tbl 348 r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 349 350 zhk = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1 351 ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb) ! proportion of bottom cell influenced by boundary layer 352 END DO 353 END DO 354 355 END SUBROUTINE sbc_isf_init 356 323 357 SUBROUTINE sbc_isf_bg03(kt) 324 !!========================================================================== 325 !! *** SUBROUTINE sbcisf_bg03 *** 326 !! add net heat and fresh water flux from ice shelf melting 327 !! into the adjacent ocean using the parameterisation by 328 !! Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean 329 !! interaction for climate models", Ocean Modelling 5(2003) 157-170. 330 !! (hereafter BG) 331 !!========================================================================== 332 !!---------------------------------------------------------------------- 333 !! sbc_isf_bg03 : routine called from sbcmod 334 !!---------------------------------------------------------------------- 335 !! 336 !! ** Purpose : Add heat and fresh water fluxes due to ice shelf melting 337 !! ** Reference : Beckmann et Goosse, 2003, Ocean Modelling 338 !! 339 !! History : 340 !! ! 06-02 (C. Wang) Original code 341 !!---------------------------------------------------------------------- 342 343 INTEGER, INTENT ( in ) :: kt 344 345 INTEGER :: ji, jj, jk !temporary integer 346 347 REAL(wp) :: zt_sum ! sum of the temperature between 200m and 600m 348 REAL(wp) :: zt_ave ! averaged temperature between 200m and 600m 349 REAL(wp) :: zt_frz ! freezing point temperature at depth z 350 REAL(wp) :: zpress ! pressure to compute the freezing point in depth 351 352 !!---------------------------------------------------------------------- 353 IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03') 354 ! 355 356 ! This test is false only in the very first time step of a run (JMM ???- Initialy build to skip 1rst year of run ) 357 DO ji = 1, jpi 358 DO jj = 1, jpj 359 jk = misfkt(ji,jj) 360 !! Initialize arrays to 0 (each step) 361 zt_sum = 0.e0_wp 362 IF ( jk .GT. 1 ) THEN 363 ! 1. -----------the average temperature between 200m and 600m --------------------- 364 DO jk = misfkt(ji,jj),misfkb(ji,jj) 365 ! freezing point temperature at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!) 366 ! after verif with UNESCO, wrong sign in BG eq. 2 367 ! Calculate freezing temperature 368 CALL eos_fzp(stbl(ji,jj), zt_frz, zpress) 369 zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) ! sum temp 370 ENDDO 371 zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value 372 373 ! 2. ------------Net heat flux and fresh water flux due to the ice shelf 374 ! For those corresponding to zonal boundary 375 qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave & 376 & / (e1t(ji,jj) * e2t(ji,jj)) * tmask(ji,jj,jk) 358 !!--------------------------------------------------------------------- 359 !! *** ROUTINE sbc_isf_bg03 *** 360 !! 361 !! ** Purpose : add net heat and fresh water flux from ice shelf melting 362 !! into the adjacent ocean 363 !! 364 !! ** Method : See reference 365 !! 366 !! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean 367 !! interaction for climate models", Ocean Modelling 5(2003) 157-170. 368 !! (hereafter BG) 369 !! History : 370 !! 06-02 (C. Wang) Original code 371 !!---------------------------------------------------------------------- 372 INTEGER, INTENT ( in ) :: kt 373 ! 374 INTEGER :: ji, jj, jk ! dummy loop index 375 INTEGER :: ik ! current level 376 REAL(wp) :: zt_sum ! sum of the temperature between 200m and 600m 377 REAL(wp) :: zt_ave ! averaged temperature between 200m and 600m 378 REAL(wp) :: zt_frz ! freezing point temperature at depth z 379 REAL(wp) :: zpress ! pressure to compute the freezing point in depth 380 !!---------------------------------------------------------------------- 381 382 IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03') 383 ! 384 DO ji = 1, jpi 385 DO jj = 1, jpj 386 ik = misfkt(ji,jj) 387 !! Initialize arrays to 0 (each step) 388 zt_sum = 0.e0_wp 389 IF ( ik > 1 ) THEN 390 ! 1. -----------the average temperature between 200m and 600m --------------------- 391 DO jk = misfkt(ji,jj),misfkb(ji,jj) 392 ! freezing point temperature at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!) 393 ! after verif with UNESCO, wrong sign in BG eq. 2 394 ! Calculate freezing temperature 395 CALL eos_fzp(stbl(ji,jj), zt_frz, zpress) 396 zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) ! sum temp 397 END DO 398 zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value 399 ! 2. ------------Net heat flux and fresh water flux due to the ice shelf 400 ! For those corresponding to zonal boundary 401 qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave & 402 & * r1_e12t(ji,jj) * tmask(ji,jj,jk) 377 403 378 fwfisf(ji,jj) = qisf(ji,jj) / lfusisf !fresh water flux kg/(m2s) 379 fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) ) 380 !add to salinity trend 381 ELSE 382 qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp 383 END IF 384 ENDDO 385 ENDDO 386 ! 387 IF( nn_timing == 1 ) CALL timing_stop('sbc_isf_bg03') 404 fwfisf(ji,jj) = qisf(ji,jj) / rlfusisf !fresh water flux kg/(m2s) 405 fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) ) 406 !add to salinity trend 407 ELSE 408 qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp 409 END IF 410 END DO 411 END DO 412 ! 413 IF( nn_timing == 1 ) CALL timing_stop('sbc_isf_bg03') 414 388 415 END SUBROUTINE sbc_isf_bg03 389 416 390 417 SUBROUTINE sbc_isf_cav( kt ) 391 418 !!--------------------------------------------------------------------- 392 419 !! *** ROUTINE sbc_isf_cav *** … … 403 430 INTEGER, INTENT(in) :: kt ! ocean time step 404 431 ! 405 REAL(wp), DIMENSION(:,:), POINTER :: zfrz 406 REAL(wp), DIMENSION(:,:), POINTER :: zgammat, zgammas 407 REAL(wp), DIMENSION(:,:), POINTER :: zfwflx, zhtflx, zhtflx_b 432 INTEGER :: ji, jj ! dummy loop indices 433 INTEGER :: nit 408 434 REAL(wp) :: zlamb1, zlamb2, zlamb3 409 435 REAL(wp) :: zeps1,zeps2,zeps3,zeps4,zeps6,zeps7 … … 411 437 REAL(wp) :: zeps = 1.e-20_wp 412 438 REAL(wp) :: zerr 413 INTEGER :: ji, jj ! dummy loop indices 414 INTEGER :: nit 439 REAL(wp), DIMENSION(:,:), POINTER :: zfrz 440 REAL(wp), DIMENSION(:,:), POINTER :: zgammat, zgammas 441 REAL(wp), DIMENSION(:,:), POINTER :: zfwflx, zhtflx, zhtflx_b 415 442 LOGICAL :: lit 416 443 !!--------------------------------------------------------------------- 417 !418 444 ! coeficient for linearisation of potential tfreez 419 445 ! Crude approximation for pressure (but commonly used) 420 zlamb1 =-0.0573_wp421 zlamb2 = 0.0832_wp422 zlamb3 =-7.53e-08* grav * rau0446 zlamb1 =-0.0573_wp 447 zlamb2 = 0.0832_wp 448 zlamb3 =-7.53e-08_wp * grav * rau0 423 449 IF( nn_timing == 1 ) CALL timing_start('sbc_isf_cav') 424 450 ! … … 427 453 428 454 ! initialisation 429 zgammat(:,:) =rn_gammat0 ; zgammas (:,:)=rn_gammas0;430 zhtflx (:,:) =0.0_wp ; zhtflx_b(:,:)=0.0_wp ;431 zfwflx (:,:) =0.0_wp455 zgammat(:,:) = rn_gammat0 ; zgammas (:,:) = rn_gammas0 456 zhtflx (:,:) = 0.0_wp ; zhtflx_b(:,:) = 0.0_wp 457 zfwflx (:,:) = 0.0_wp 432 458 433 459 ! compute ice shelf melting 434 460 nit = 1 ; lit = .TRUE. 435 461 DO WHILE ( lit ) ! maybe just a constant number of iteration as in blk_core is fine 436 IF (nn_isfblk == 1) THEN437 !ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006)462 SELECT CASE ( nn_isfblk ) 463 CASE ( 1 ) ! ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006) 438 464 ! Calculate freezing temperature 439 465 CALL eos_fzp( stbl(:,:), zfrz(:,:), risfdep(:,:) ) … … 446 472 DO ji = 1, jpi 447 473 zhtflx(ji,jj) = zgammat(ji,jj)*rcp*rau0*(ttbl(ji,jj)-zfrz(ji,jj)) 448 zfwflx(ji,jj) = - zhtflx(ji,jj)/ lfusisf474 zfwflx(ji,jj) = - zhtflx(ji,jj)/rlfusisf 449 475 END DO 450 476 END DO … … 454 480 fwfisf(:,:) = zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 455 481 456 ELSE IF (nn_isfblk == 2 ) THEN 457 ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015) 482 CASE ( 2 ) ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015) 458 483 ! compute gammat every where (2d) 459 484 CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) … … 464 489 DO ji = 1, jpi 465 490 ! compute coeficient to solve the 2nd order equation 466 zeps1 =rcp*rau0*zgammat(ji,jj)467 zeps2 =lfusisf*rau0*zgammas(ji,jj)468 zeps3 =rhoisf*rcpi*kappa/MAX(risfdep(ji,jj),zeps)469 zeps4 =zlamb2+zlamb3*risfdep(ji,jj)470 zeps6 =zeps4-ttbl(ji,jj)471 zeps7 =zeps4-tsurf472 zaqe =zlamb1 * (zeps1 + zeps3)473 zaqer =0.5/MIN(zaqe,-zeps)474 zbqe =zeps1*zeps6+zeps3*zeps7-zeps2475 zcqe =zeps2*stbl(ji,jj)476 zdis =zbqe*zbqe-4.0*zaqe*zcqe491 zeps1 = rcp*rau0*zgammat(ji,jj) 492 zeps2 = rlfusisf*rau0*zgammas(ji,jj) 493 zeps3 = rhoisf*rcpi*rkappa/MAX(risfdep(ji,jj),zeps) 494 zeps4 = zlamb2+zlamb3*risfdep(ji,jj) 495 zeps6 = zeps4-ttbl(ji,jj) 496 zeps7 = zeps4-tsurf 497 zaqe = zlamb1 * (zeps1 + zeps3) 498 zaqer = 0.5_wp/MIN(zaqe,-zeps) 499 zbqe = zeps1*zeps6+zeps3*zeps7-zeps2 500 zcqe = zeps2*stbl(ji,jj) 501 zdis = zbqe*zbqe-4.0_wp*zaqe*zcqe 477 502 478 503 ! Presumably zdis can never be negative because gammas is very small compared to gammat 479 504 ! compute s freeze 480 505 zsfrz=(-zbqe-SQRT(zdis))*zaqer 481 IF ( zsfrz .LT.0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer506 IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer 482 507 483 508 ! compute t freeze (eq. 22) … … 496 521 fwfisf(:,:) = zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 497 522 498 END IF523 END SELECT 499 524 500 525 ! define if we need to iterate (nn_gammablk 0/1 do not need iteration) 501 IF ( nn_gammablk .LT.2 ) THEN ; lit = .FALSE.526 IF ( nn_gammablk < 2 ) THEN ; lit = .FALSE. 502 527 ELSE 503 528 ! check total number of iteration 504 IF (nit .GE.100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )505 ELSE 506 END IF529 IF (nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 530 ELSE ; nit = nit + 1 531 END IF 507 532 508 533 ! compute error between 2 iterations 509 534 ! if needed save gammat and compute zhtflx_b for next iteration 510 535 zerr = MAXVAL(ABS(zhtflx-zhtflx_b)) 511 IF ( zerr .LE. 0.01) THEN ; lit = .FALSE.512 ELSE ; zhtflx_b(:,:) = zhtflx(:,:)513 END IF536 IF ( zerr <= 0.01_wp ) THEN ; lit = .FALSE. 537 ELSE ; zhtflx_b(:,:) = zhtflx(:,:) 538 END IF 514 539 END IF 515 540 END DO 516 541 ! 517 ! output 518 IF( iom_use('isfgammat') ) CALL iom_put('isfgammat', zgammat) 519 IF( iom_use('isfgammas') ) CALL iom_put('isfgammas', zgammas) 542 ! output see JMM comment above 543 ! IF( iom_use('isfgammat') ) CALL iom_put('isfgammat', zgammat) 544 ! IF( iom_use('isfgammas') ) CALL iom_put('isfgammas', zgammas) 545 CALL iom_put('isfgammat', zgammat) 546 CALL iom_put('isfgammas', zgammas) 520 547 ! 521 548 CALL wrk_dealloc( jpi,jpj, zfrz , zgammat, zgammas ) … … 526 553 END SUBROUTINE sbc_isf_cav 527 554 528 SUBROUTINE sbc_isf_gammats( gt, gs, zqhisf, zqwisf )555 SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf ) 529 556 !!---------------------------------------------------------------------- 530 557 !! ** Purpose : compute the coefficient echange for heat flux … … 535 562 !! Jenkins et al., 2010, JPO, p2298-2312 536 563 !!--------------------------------------------------------------------- 537 REAL(wp), DIMENSION(:,:), INTENT(out) :: gt,gs538 REAL(wp), DIMENSION(:,:), INTENT(in ) :: zqhisf, zqwisf539 564 REAL(wp), DIMENSION(:,:), INTENT(out) :: pgt, pgs 565 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pqhisf, pqwisf 566 ! 540 567 INTEGER :: ikt 541 INTEGER :: ji, jj! loop index568 INTEGER :: ji, jj ! loop index 542 569 REAL(wp), DIMENSION(:,:), POINTER :: zustar ! U, V at T point and friction velocity 543 570 REAL(wp) :: zdku, zdkv ! U, V shear … … 551 578 REAL(wp) :: zdep 552 579 REAL(wp) :: zeps = 1.0e-20_wp 553 REAL(wp), PARAMETER :: zxsiN = 0.052 ! dimensionless constant554 REAL(wp), PARAMETER :: znu = 1.95e-6 ! kinamatic viscosity of sea water (m2.s-1)580 REAL(wp), PARAMETER :: zxsiN = 0.052_wp ! dimensionless constant 581 REAL(wp), PARAMETER :: znu = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1) 555 582 REAL(wp), DIMENSION(2) :: zts, zab 556 583 !!--------------------------------------------------------------------- 557 584 CALL wrk_alloc( jpi,jpj, zustar ) 558 585 ! 559 IF ( nn_gammablk == 0 ) THEN 560 !! gamma is constant (specified in namelist) 561 !! ISOMIP formulation (Hunter et al, 2006) 562 gt(:,:) = rn_gammat0 563 gs(:,:) = rn_gammas0 564 565 ELSE IF ( nn_gammablk == 1 ) THEN 566 !! gamma is assume to be proportional to u* 567 !! Jenkins et al., 2010, JPO, p2298-2312 568 !! Adopted by Asay-Davis et al. (2015) 569 570 !! compute ustar (eq. 24) 586 SELECT CASE ( nn_gammablk ) 587 CASE ( 0 ) ! gamma is constant (specified in namelist) 588 !! ISOMIP formulation (Hunter et al, 2006) 589 pgt(:,:) = rn_gammat0 590 pgs(:,:) = rn_gammas0 591 592 CASE ( 1 ) ! gamma is assume to be proportional to u* 593 !! Jenkins et al., 2010, JPO, p2298-2312 594 !! Adopted by Asay-Davis et al. (2015) 595 596 !! compute ustar (eq. 24) 571 597 zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) ) 572 598 573 !! Compute gammats574 gt(:,:) = zustar(:,:) * rn_gammat0575 gs(:,:) = zustar(:,:) * rn_gammas0599 !! Compute gammats 600 pgt(:,:) = zustar(:,:) * rn_gammat0 601 pgs(:,:) = zustar(:,:) * rn_gammas0 576 602 577 ELSE IF ( nn_gammablk == 2 ) THEN 578 !! gamma depends of stability of boundary layer 579 !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14 580 !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO) 581 !! compute ustar 603 CASE ( 2 ) ! gamma depends of stability of boundary layer 604 !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14 605 !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO) 606 !! compute ustar 582 607 zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) ) 583 608 584 !! compute Pr and Sc number (can be improved)585 zPr = 13.8 586 zSc = 2432.0 587 588 !! compute gamma mole589 zgmolet = 12.5 * zPr ** (2.0/3.0) - 6.0590 zgmoles = 12.5 * zSc ** (2.0/3.0) -6.0591 592 !! compute gamma609 !! compute Pr and Sc number (can be improved) 610 zPr = 13.8_wp 611 zSc = 2432.0_wp 612 613 !! compute gamma mole 614 zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp 615 zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp 616 617 !! compute gamma 593 618 DO ji=2,jpi 594 619 DO jj=2,jpj … … 596 621 597 622 IF (zustar(ji,jj) == 0._wp) THEN ! only for kt = 1 I think 598 gt = rn_gammat0599 gs = rn_gammas0623 pgt = rn_gammat0 624 pgs = rn_gammas0 600 625 ELSE 601 !! compute Rc number (as done in zdfric.F90)602 zcoef = 0.5 / fse3w(ji,jj,ikt)626 !! compute Rc number (as done in zdfric.F90) 627 zcoef = 0.5_wp / fse3w(ji,jj,ikt) 603 628 ! ! shear of horizontal velocity 604 629 zdku = zcoef * ( un(ji-1,jj ,ikt ) + un(ji,jj,ikt ) & … … 609 634 zRc = rn2(ji,jj,ikt+1) / MAX( zdku*zdku + zdkv*zdkv, zeps ) 610 635 611 !! compute bouyancy636 !! compute bouyancy 612 637 zts(jp_tem) = ttbl(ji,jj) 613 638 zts(jp_sal) = stbl(ji,jj) … … 616 641 CALL eos_rab( zts, zdep, zab ) 617 642 ! 618 !! compute length scale619 zbuofdep = grav * ( zab(jp_tem) * zqhisf(ji,jj) - zab(jp_sal) * zqwisf(ji,jj) ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!620 621 !! compute Monin Obukov Length643 !! compute length scale 644 zbuofdep = grav * ( zab(jp_tem) * pqhisf(ji,jj) - zab(jp_sal) * pqwisf(ji,jj) ) !!!!!!!!!!!!!!!!!!!!!!!!!!!! 645 646 !! compute Monin Obukov Length 622 647 ! Maximum boundary layer depth 623 648 zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) - 0.001_wp … … 626 651 zmols = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt) 627 652 628 !! compute eta* (stability parameter)653 !! compute eta* (stability parameter) 629 654 zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0_wp))) 630 655 631 !! compute the sublayer thickness656 !! compute the sublayer thickness 632 657 zhnu = 5 * znu / zustar(ji,jj) 633 658 634 !! compute gamma turb659 !! compute gamma turb 635 660 zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) & 636 661 & + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn 637 662 638 !! compute gammats639 gt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet)640 gs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles)663 !! compute gammats 664 pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet) 665 pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles) 641 666 END IF 642 667 END DO 643 668 END DO 644 CALL lbc_lnk( gt(:,:),'T',1.)645 CALL lbc_lnk( gs(:,:),'T',1.)646 END IF669 CALL lbc_lnk(pgt(:,:),'T',1.) 670 CALL lbc_lnk(pgs(:,:),'T',1.) 671 END SELECT 647 672 CALL wrk_dealloc( jpi,jpj, zustar ) 648 673 649 END SUBROUTINE 650 651 SUBROUTINE sbc_isf_tbl( varin, varout, cptin )674 END SUBROUTINE sbc_isf_gammats 675 676 SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin ) 652 677 !!---------------------------------------------------------------------- 653 678 !! *** SUBROUTINE sbc_isf_tbl *** … … 656 681 !! 657 682 !!---------------------------------------------------------------------- 658 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: varin 659 REAL(wp), DIMENSION(:,:) , INTENT(out):: varout 660 661 CHARACTER(len=1), INTENT(in) :: cptin ! point of variable in/out 662 683 REAL(wp), DIMENSION(:,:,:), INTENT( in ) :: pvarin 684 REAL(wp), DIMENSION(:,:) , INTENT( out ) :: pvarout 685 CHARACTER(len=1), INTENT( in ) :: cd_ptin ! point of variable in/out 686 ! 663 687 REAL(wp) :: ze3, zhk 664 688 REAL(wp), DIMENSION(:,:), POINTER :: zhisf_tbl ! thickness of the tbl 665 689 666 INTEGER :: ji, jj,jk! loop index667 INTEGER :: ikt, ikb ! top and bottom index of the tbl668 690 INTEGER :: ji, jj, jk ! loop index 691 INTEGER :: ikt, ikb ! top and bottom index of the tbl 692 !!---------------------------------------------------------------------- 669 693 ! allocation 670 694 CALL wrk_alloc( jpi,jpj, zhisf_tbl) 671 695 672 696 ! initialisation 673 varout(:,:)=0._wp697 pvarout(:,:)=0._wp 674 698 675 ! compute U in the top boundary layer at T- point676 IF (cptin == 'U') THEN699 SELECT CASE ( cd_ptin ) 700 CASE ( 'U' ) ! compute U in the top boundary layer at T- point 677 701 DO jj = 1,jpj 678 702 DO ji = 1,jpi 679 703 ikt = miku(ji,jj) ; ikb = miku(ji,jj) 680 ! thickness of boundary layer at least the top level thickness704 ! thickness of boundary layer at least the top level thickness 681 705 zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3u_n(ji,jj,ikt)) 682 706 683 ! determine the deepest level influenced by the boundary layer707 ! determine the deepest level influenced by the boundary layer 684 708 DO jk = ikt+1, mbku(ji,jj) 685 IF ( (SUM(fse3u_n(ji,jj,ikt:jk-1)) .LT.zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk709 IF ( (SUM(fse3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk 686 710 END DO 687 711 zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3u_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 688 712 689 ! level fully include in the ice shelf boundary layer713 ! level fully include in the ice shelf boundary layer 690 714 DO jk = ikt, ikb - 1 691 715 ze3 = fse3u_n(ji,jj,jk) 692 varout(ji,jj) = varout(ji,jj) +varin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3693 END DO 694 695 ! level partially include in ice shelf boundary layer716 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 717 END DO 718 719 ! level partially include in ice shelf boundary layer 696 720 zhk = SUM( fse3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 697 varout(ji,jj) = varout(ji,jj) +varin(ji,jj,ikb) * (1._wp - zhk)721 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 698 722 END DO 699 723 END DO 700 724 DO jj = 2,jpj 701 725 DO ji = 2,jpi 702 varout(ji,jj) = 0.5_wp * (varout(ji,jj) + varout(ji-1,jj)) 703 END DO 704 END DO 705 CALL lbc_lnk(varout,'T',-1.) 706 END IF 707 708 ! compute V in the top boundary layer at T- point 709 IF (cptin == 'V') THEN 726 pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji-1,jj)) 727 END DO 728 END DO 729 CALL lbc_lnk(pvarout,'T',-1.) 730 731 CASE ( 'V' ) ! compute V in the top boundary layer at T- point 710 732 DO jj = 1,jpj 711 733 DO ji = 1,jpi 712 734 ikt = mikv(ji,jj) ; ikb = mikv(ji,jj) 713 ! thickness of boundary layer at least the top level thickness735 ! thickness of boundary layer at least the top level thickness 714 736 zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3v_n(ji,jj,ikt)) 715 737 716 ! determine the deepest level influenced by the boundary layer738 ! determine the deepest level influenced by the boundary layer 717 739 DO jk = ikt+1, mbkv(ji,jj) 718 IF ( (SUM(fse3v_n(ji,jj,ikt:jk-1)) .LT.zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk740 IF ( (SUM(fse3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk 719 741 END DO 720 742 zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3v_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 721 743 722 ! level fully include in the ice shelf boundary layer744 ! level fully include in the ice shelf boundary layer 723 745 DO jk = ikt, ikb - 1 724 746 ze3 = fse3v_n(ji,jj,jk) 725 varout(ji,jj) = varout(ji,jj) +varin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3726 END DO 727 728 ! level partially include in ice shelf boundary layer747 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 748 END DO 749 750 ! level partially include in ice shelf boundary layer 729 751 zhk = SUM( fse3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 730 varout(ji,jj) = varout(ji,jj) +varin(ji,jj,ikb) * (1._wp - zhk)752 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 731 753 END DO 732 754 END DO 733 755 DO jj = 2,jpj 734 756 DO ji = 2,jpi 735 varout(ji,jj) = 0.5_wp * (varout(ji,jj) + varout(ji,jj-1)) 736 END DO 737 END DO 738 CALL lbc_lnk(varout,'T',-1.) 739 END IF 740 741 ! compute T in the top boundary layer at T- point 742 IF (cptin == 'T') THEN 757 pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji,jj-1)) 758 END DO 759 END DO 760 CALL lbc_lnk(pvarout,'T',-1.) 761 762 CASE ( 'T' ) ! compute T in the top boundary layer at T- point 743 763 DO jj = 1,jpj 744 764 DO ji = 1,jpi … … 746 766 ikb = misfkb(ji,jj) 747 767 748 ! level fully include in the ice shelf boundary layer768 ! level fully include in the ice shelf boundary layer 749 769 DO jk = ikt, ikb - 1 750 770 ze3 = fse3t_n(ji,jj,jk) 751 varout(ji,jj) = varout(ji,jj) +varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3752 END DO 753 754 ! level partially include in ice shelf boundary layer771 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3 772 END DO 773 774 ! level partially include in ice shelf boundary layer 755 775 zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) 756 varout(ji,jj) = varout(ji,jj) +varin(ji,jj,ikb) * (1._wp - zhk)757 END DO 758 END DO 759 END IF776 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 777 END DO 778 END DO 779 END SELECT 760 780 761 781 ! mask mean tbl value 762 varout(:,:) =varout(:,:) * ssmask(:,:)782 pvarout(:,:) = pvarout(:,:) * ssmask(:,:) 763 783 764 784 ! deallocation … … 780 800 !! ** Action : phdivn decreased by the runoff inflow 781 801 !!---------------------------------------------------------------------- 782 REAL(wp), DIMENSION(:,:,:), INTENT( inout) :: phdivn ! horizontal divergence783 ! !802 REAL(wp), DIMENSION(:,:,:), INTENT( inout ) :: phdivn ! horizontal divergence 803 ! 784 804 INTEGER :: ji, jj, jk ! dummy loop indices 785 805 INTEGER :: ikt, ikb … … 790 810 zfact = 0.5_wp 791 811 ! 792 IF ( lk_vvl) THEN ! need to re compute level distribution of isf fresh water812 IF ( lk_vvl ) THEN ! need to re compute level distribution of isf fresh water 793 813 DO jj = 1,jpj 794 814 DO ji = 1,jpi … … 800 820 ! determine the deepest level influenced by the boundary layer 801 821 DO jk = ikt+1, mbkt(ji,jj) 802 IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT.rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk822 IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 803 823 END DO 804 824 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. … … 820 840 DO jk = ikt, ikb - 1 821 841 phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) & 822 & 842 & * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact 823 843 END DO 824 844 ! level partially include in ice shelf boundary layer 825 845 phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) & 826 &+ fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj)846 & + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 827 847 END DO 828 848 END DO -
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
r5621 r5944 141 141 IF( ln_rnf_tem ) THEN ! use runoffs temperature data 142 142 rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0 143 CALL eos_fzp( sss_m(:,:), ztfrz(:,:) ) 143 144 WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999._wp ) ! if missing data value use SST as runoffs temperature 144 145 rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0 145 146 END WHERE 146 147 WHERE( sf_t_rnf(1)%fnow(:,:,1) == -222._wp ) ! where fwf comes from melting of ice shelves or iceberg 147 ztfrz(:,:) = -1.9 !tfreez( sss_m(:,:) ) !PM to be discuss (trouble if sensitivity study) 148 rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * lfusisf * r1_rau0_rcp 148 rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * rlfusisf * r1_rau0_rcp 149 149 END WHERE 150 150 ELSE ! use SST as runoffs temperature
Note: See TracChangeset
for help on using the changeset viewer.