Changeset 5066 for branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO
- Timestamp:
- 2015-02-06T17:02:20+01:00 (9 years ago)
- Location:
- branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 1 added
- 9 edited
-
DOM/dom_oce.F90 (modified) (4 diffs)
-
DOM/domain.F90 (modified) (1 diff)
-
DOM/domvvl.F90 (modified) (3 diffs)
-
DOM/domzgr.F90 (modified) (17 diffs)
-
DYN/dynhpg.F90 (modified) (20 diffs)
-
DYN/dynspg_ts.F90 (modified) (17 diffs)
-
DYN/sshwzv.F90 (modified) (4 diffs)
-
DYN/wadlmt.F90 (added)
-
nemogcm.F90 (modified) (5 diffs)
-
step_oce.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r4488 r5066 10 10 !! 3.5 ! 2012 (S. Mocavero, I. Epicoco) Add arrays associated 11 11 !! to the optimization of BDY communications 12 !! 3.6.?! 2014 (H. Liu) Add arrays associated with Wetting and Drying 12 13 !!---------------------------------------------------------------------- 13 14 … … 262 263 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: nicoa, njcoa !: ??? 263 264 #endif 265 266 !!---------------------------------------------------------------------- 267 !! critical depths,limiters,and masks for Wetting and Drying 268 !! --------------------------------------------------------------------- 269 270 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wdmask !: u- and v- limiter 271 272 LOGICAL, PUBLIC, SAVE :: ln_wd !: key to turn on/off wetting/drying (T: on, F: off) 273 REAL(wp), PUBLIC, SAVE :: rn_wdmin1 !: minimum water depth on dried cells 274 REAL(wp), PUBLIC, SAVE :: rn_wdmin2 !: tolerrance of minimum water depth on dried cells 275 REAL(wp), PUBLIC, SAVE :: rn_wdld !: land elevation below which wetting/drying will be considered 276 INTEGER , PUBLIC, SAVE :: nn_wdit !: maximum number of iteration for W/D limiter 264 277 265 278 !!---------------------------------------------------------------------- … … 325 338 INTEGER FUNCTION dom_oce_alloc() 326 339 !!---------------------------------------------------------------------- 327 INTEGER, DIMENSION(1 1) :: ierr340 INTEGER, DIMENSION(12) :: ierr 328 341 !!---------------------------------------------------------------------- 329 342 ierr(:) = 0 … … 390 403 #endif 391 404 ! 405 IF(ln_wd) & 406 ALLOCATE( wdmask(jpi,jpj), STAT=ierr(12) ) 407 ! 392 408 dom_oce_alloc = MAXVAL(ierr) 393 409 ! -
branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
r4624 r5066 86 86 CALL dom_zgr ! Vertical mesh and bathymetry 87 87 CALL dom_msk ! Masks 88 IF( ln_sco ) CALL dom_stiff ! Maximum stiffness ratio/hydrostatic consistency 88 IF( ln_sco.AND.(.NOT.ln_wd) ) & 89 & CALL dom_stiff ! Maximum stiffness ratio/hydrostatic consistency 89 90 ! 90 91 ht_0(:,:) = 0.0_wp ! Reference ocean depth at T-points -
branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r4624 r5066 145 145 fse3t_a(:,:,jpk) = e3t_0(:,:,jpk) 146 146 147 IF(ln_wd) THEN 148 DO jj = 1, jpj 149 DO ji = 1, jpi 150 IF(mbathy(ji,jj) == 2 .AND. e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1) THEN 151 fse3t_a(ji,jj,1:2) = 0.5_wp * rn_wdmin1 152 END IF 153 ENDDO 154 ENDDO 155 END IF 156 147 157 ! Reconstruction of all vertical scale factors at now and before time steps 148 158 ! ============================================================================= … … 793 803 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 794 804 !! * Local declarations 795 INTEGER :: j k805 INTEGER :: ji, jj, jk 796 806 INTEGER :: id1, id2, id3, id4, id5 ! local integers 797 807 !!---------------------------------------------------------------------- … … 867 877 fse3t_n(:,:,:) = e3t_0(:,:,:) 868 878 sshn(:,:) = 0.0_wp 879 880 IF(ln_wd) THEN 881 DO jj = 1, jpj 882 DO ji = 1, jpi 883 !IF(e3t_0(ji,jj,1) < 0._wp) THEN 884 IF(mbathy(ji,jj) == 2 .AND. e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1) THEN 885 fse3t_b(ji,jj,1:2) = 0.5_wp * rn_wdmin1 886 fse3t_n(ji,jj,1:2) = 0.5_wp * rn_wdmin1 887 sshb(ji,jj) = rn_wdmin1 - bathy(ji,jj) 888 sshn(ji,jj) = rn_wdmin1 - bathy(ji,jj) 889 ENDIF 890 ENDDO 891 ENDDO 892 END IF 893 869 894 IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN 870 895 tilde_e3t_b(:,:,:) = 0.0_wp -
branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r4754 r5066 17 17 !! 3.4 ! 2012-08 (J. Siddorn) added Siddorn and Furner stretching function 18 18 !! 3.4 ! 2012-12 (R. Bourdalle-Badie and G. Reffray) modify C1D case 19 !! 3.6 ! 2014-09 (H. Liu) Modifications for Wetting/Drying 19 20 !!---------------------------------------------------------------------- 20 21 … … 1121 1122 REAL(wp) :: zrmax, ztaper ! temporary scalars 1122 1123 REAL(wp) :: zrfact 1124 REAL(wp) :: zsmth 1123 1125 ! 1124 1126 REAL(wp), POINTER, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 … … 1176 1178 bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) ) 1177 1179 1178 DO jj = 1, jpj 1179 DO ji = 1, jpi 1180 IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) 1181 END DO 1182 END DO 1180 IF (.NOT. ln_wd) THEN 1181 DO jj = 1, jpj 1182 DO ji = 1, jpi 1183 IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) 1184 END DO 1185 END DO 1186 ENDIF 1183 1187 ! ! ============================= 1184 1188 ! ! Define the envelop bathymetry (hbatt) … … 1187 1191 zenv(:,:) = bathy(:,:) 1188 1192 ! 1189 ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing 1190 DO jj = 1, jpj 1191 DO ji = 1, jpi 1192 IF( bathy(ji,jj) == 0._wp ) THEN 1193 iip1 = MIN( ji+1, jpi ) 1194 ijp1 = MIN( jj+1, jpj ) 1195 iim1 = MAX( ji-1, 1 ) 1196 ijm1 = MAX( jj-1, 1 ) 1197 IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) + & 1198 & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 1199 zenv(ji,jj) = rn_sbot_min 1193 IF (.NOT. ln_wd) THEN 1194 ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing 1195 DO jj = 1, jpj 1196 DO ji = 1, jpi 1197 IF( bathy(ji,jj) == 0._wp ) THEN 1198 iip1 = MIN( ji+1, jpi ) 1199 ijp1 = MIN( jj+1, jpj ) 1200 iim1 = MAX( ji-1, 1 ) 1201 ijm1 = MAX( jj-1, 1 ) 1202 IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) + & 1203 & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 1204 zenv(ji,jj) = rn_sbot_min 1205 ENDIF 1200 1206 ENDIF 1201 END IF1202 1203 END DO1207 END DO 1208 END DO 1209 END IF 1204 1210 ! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero 1205 1211 CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' ) … … 1241 1247 zri(:,:) = 0._wp 1242 1248 zrj(:,:) = 0._wp 1249 zsmth = 0._wp 1243 1250 DO jj = 1, nlcj 1244 1251 DO ji = 1, nlci 1245 1252 iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi) 1246 1253 ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj) 1247 IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN1254 IF( (zenv(ji,jj) > zsmth) .AND. (zenv(iip1,jj) > zsmth)) THEN 1248 1255 zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) 1249 1256 END IF 1250 IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN1257 IF( (zenv(ji,jj) > zsmth) .AND. (zenv(ji,ijp1) > zsmth)) THEN 1251 1258 zrj(ji,jj) = ( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) 1252 1259 END IF … … 1271 1278 END DO ! End loop ! 1272 1279 ! ! ================ ! 1273 DO jj = 1, jpj 1274 DO ji = 1, jpi 1275 zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings 1276 END DO 1277 END DO 1280 IF(ln_wd) THEN 1281 DO jj = 1, jpj 1282 DO ji = 1, jpi 1283 zenv(ji,jj) = MAX( zenv(ji,jj), -rn_wdld ) ! fill out land bathy data 1284 END DO 1285 END DO 1286 ELSE 1287 DO jj = 1, jpj 1288 DO ji = 1, jpi 1289 zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings 1290 END DO 1291 END DO 1292 END IF 1278 1293 ! 1279 1294 ! Envelope bathymetry saved in hbatt … … 1305 1320 IF(lwp) THEN 1306 1321 WRITE(numout,*) 1307 WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 1322 IF (.NOT. ln_wd ) THEN 1323 WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 1324 ELSE 1325 WRITE(numout,*) ' zgr_sco: minimum positive depth of the envelop topography set to : ', rn_sbot_min 1326 WRITE(numout,*) ' zgr_sco: minimum negative depth of the envelop topography set to : ', -rn_wdld 1327 ENDIF 1308 1328 ENDIF 1309 1329 hbatu(:,:) = rn_sbot_min … … 1318 1338 END DO 1319 1339 END DO 1340 1341 IF(ln_wd) THEN !avoid the zero depth on T- (U-,V-,F-) points 1342 DO jj = 1, jpj 1343 DO ji = 1, jpi 1344 IF(ABS(hbatt(ji,jj)) < rn_wdmin1) & 1345 & hbatt(ji,jj) = SIGN(1._wp, hbatt(ji,jj)) * rn_wdmin1 1346 IF(ABS(hbatu(ji,jj)) < rn_wdmin1) & 1347 & hbatu(ji,jj) = SIGN(1._wp, hbatu(ji,jj)) * rn_wdmin1 1348 IF(ABS(hbatv(ji,jj)) < rn_wdmin1) & 1349 & hbatv(ji,jj) = SIGN(1._wp, hbatv(ji,jj)) * rn_wdmin1 1350 IF(ABS(hbatf(ji,jj)) < rn_wdmin1) & 1351 & hbatf(ji,jj) = SIGN(1._wp, hbatf(ji,jj)) * rn_wdmin1 1352 END DO 1353 END DO 1354 END IF 1320 1355 ! 1321 1356 ! Apply lateral boundary condition … … 1325 1360 DO ji = 1, jpi 1326 1361 IF( hbatu(ji,jj) == 0._wp ) THEN 1362 !No worries about the following line when ln_wd == .true. 1327 1363 IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min 1328 1364 IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj) … … 1350 1386 1351 1387 !!bug: key_helsinki a verifer 1352 hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 1353 hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 1354 hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 1355 hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 1388 1389 IF (.NOT. ln_wd) THEN 1390 hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 1391 hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 1392 hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 1393 hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 1394 END IF 1356 1395 1357 1396 IF( nprint == 1 .AND. lwp ) THEN … … 1395 1434 CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 1396 1435 1397 fsdepw(:,:,:) = gdepw_0 (:,:,:) 1398 fsde3w(:,:,:) = gdep3w_0(:,:,:) 1399 ! 1400 where (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.0 1401 where (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.0 1402 where (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.0 1403 where (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.0 1404 where (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.0 1405 where (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.0 1406 where (e3vw_0 (:,:,:).eq.0.0) e3vw_0(:,:,:) = 1.0 1436 !fsdepw(:,:,:) = gdepw_0 (:,:,:) 1437 !fsde3w(:,:,:) = gdep3w_0(:,:,:) 1438 ! 1439 IF (.NOT. ln_wd) THEN 1440 where (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.0 1441 where (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.0 1442 where (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.0 1443 where (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.0 1444 where (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.0 1445 where (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.0 1446 where (e3vw_0 (:,:,:).eq.0.0) e3vw_0(:,:,:) = 1.0 1447 END IF 1407 1448 1408 1449 #if defined key_agrif … … 1449 1490 IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 1450 1491 #endif 1451 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 1452 END DO 1492 END DO 1493 IF(ln_wd) THEN 1494 IF( scobot(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN 1495 mbathy(ji,jj) = 0 1496 ELSEIF(scobot(ji,jj) <= rn_wdmin1) THEN 1497 #if key_surge 1498 mbathy(ji,jj) = 1 1499 #else 1500 mbathy(ji,jj) = 2 1501 #endif 1502 ENDIF 1503 ELSE 1504 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 1505 END IF 1453 1506 END DO 1454 1507 END DO … … 1508 1561 DO jj = 1, jpj 1509 1562 1510 IF( hbatt(ji,jj) > 0._wp) THEN1563 IF( bathy(ji,jj) > 0._wp) THEN 1511 1564 DO jk = 1, mbathy(ji,jj) 1512 1565 ! check coordinate is monotonically increasing 1513 IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 1514 WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 1515 WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 1566 !RF test... 1567 !IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 1568 ! WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 1569 ! WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 1570 IF (fse3w(ji,jj,jk) < 0._wp .OR. fse3t(ji,jj,jk) < 0._wp ) THEN 1571 WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t < 0 at point (i,j,k)= ', ji, jj, jk 1572 WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t < 0 at point (i,j,k)= ', ji, jj, jk 1516 1573 WRITE(numout,*) 'e3w',fse3w(ji,jj,:) 1517 1574 WRITE(numout,*) 'e3t',fse3t(ji,jj,:) … … 1572 1629 INTEGER :: ji, jj, jk ! dummy loop argument 1573 1630 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 1631 REAL(wp) :: ztmpu, ztmpv, ztmpf 1632 REAL(wp) :: ztmpu1, ztmpv1, ztmpf1 1574 1633 ! 1575 1634 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 … … 1625 1684 DO ji = 1, jpim1 1626 1685 DO jj = 1, jpjm1 1686 ! extended for Wetting/Drying case 1687 ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj) 1688 ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1) 1689 ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) 1690 ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj) 1691 ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1) 1692 ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * & 1693 & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 1627 1694 DO jk = 1, jpk 1628 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 1629 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1630 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 1631 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1632 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) & 1633 & + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 1634 & / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 1635 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 1636 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1637 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 1638 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1695 IF((ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1).AND.ln_wd) THEN 1696 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 1697 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 1698 ELSE 1699 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 1700 & / ztmpu 1701 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 1702 & / ztmpu 1703 END IF 1704 1705 IF((ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1).AND.ln_wd) THEN 1706 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 1707 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 1708 ELSE 1709 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 1710 & / ztmpv 1711 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 1712 & / ztmpv 1713 END IF 1714 1715 IF((ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1).AND.ln_wd) THEN 1716 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) & 1717 & + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 1718 ELSE 1719 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) & 1720 & + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 1721 & / ztmpf 1722 END IF 1723 1639 1724 ! 1640 1725 e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) … … 1676 1761 REAL(wp) :: zsmth ! smoothing around critical depth 1677 1762 REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space 1763 REAL(wp) :: ztmpu, ztmpv, ztmpf 1764 REAL(wp) :: ztmpu1, ztmpv1, ztmpf1 1678 1765 ! 1679 1766 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 … … 1754 1841 DO jj=1,jpj-1 1755 1842 1756 DO jk = 1, jpk 1757 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / & 1758 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1759 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / & 1760 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1761 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) + & 1762 hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / & 1763 ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 1764 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / & 1765 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1766 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / & 1767 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1843 ! extend to suit for Wetting/Drying case 1844 ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj) 1845 ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1) 1846 ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) 1847 ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj) 1848 ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1) 1849 ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * & 1850 & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 1851 DO jk = 1, jpk 1852 IF((ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1).AND.ln_wd) THEN 1853 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 1854 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 1855 ELSE 1856 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 1857 & / ztmpu 1858 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 1859 & / ztmpu 1860 END IF 1861 1862 IF((ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1).AND.ln_wd) THEN 1863 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 1864 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 1865 ELSE 1866 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 1867 & / ztmpv 1868 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 1869 & / ztmpv 1870 END IF 1871 1872 IF((ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1).AND.ln_wd) THEN 1873 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) & 1874 & + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 1875 ELSE 1876 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) & 1877 & + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 1878 & / ztmpf 1879 END IF 1768 1880 1769 1881 e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk) -
branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r4794 r5066 16 16 !! 3.4 ! 2011-11 (H. Liu) hpg_prj: Original code for s-coordinates 17 17 !! ! (A. Coward) suppression of hel, wdj and rot options 18 !! 3.6? ! 2014-09 (H. Liu) add Wetting/Drying pressure filter 18 19 !!---------------------------------------------------------------------- 19 20 … … 369 370 !! 370 371 INTEGER :: ji, jj, jk ! dummy loop indices 371 REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars 372 REAL(wp) :: zcoef0, zuap, zvap, znad, ztmp ! temporary scalars 373 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 372 374 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 375 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter 373 376 !!---------------------------------------------------------------------- 374 377 ! 375 378 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 379 IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 376 380 ! 377 381 IF( kt == nit000 ) THEN … … 386 390 IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume 387 391 ELSE ; znad = 0._wp ! Fixed volume 392 ENDIF 393 394 IF(ln_wd) THEN 395 DO jj = 2, jpjm1 396 DO ji = 2, jpim1 397 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 398 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 399 & rn_wdmin1 + rn_wdmin2 400 401 IF(ll_tmp1) THEN 402 zcpx(ji,jj) = 1.0_wp 403 ELSE IF(ll_tmp2) THEN 404 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 405 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 406 & (sshn(ji+1,jj) - sshn(ji,jj))) 407 ELSE 408 zcpx(ji,jj) = 0._wp 409 END IF 410 411 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 412 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 413 & rn_wdmin1 + rn_wdmin2 414 415 IF(ll_tmp1) THEN 416 zcpy(ji,jj) = 1.0_wp 417 ELSE IF(ll_tmp2) THEN 418 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 419 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 420 & (sshn(ji,jj+1) - sshn(ji,jj))) 421 ELSE 422 zcpy(ji,jj) = 0._wp 423 END IF 424 END DO 425 END DO 426 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 388 427 ENDIF 389 428 … … 402 441 zvap = -zcoef0 * ( 2._wp * znad ) & 403 442 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 443 444 IF(ln_wd) THEN 445 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 446 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 447 zuap = zuap * zcpx(ji,jj) 448 zvap = zvap * zcpy(ji,jj) 449 ENDIF 450 404 451 ! add to the general momentum trend 405 452 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap … … 424 471 zvap = -zcoef0 * ( 2._wp * znad ) & 425 472 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 473 474 IF(ln_wd) THEN 475 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 476 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 477 zuap = zuap * zcpx(ji,jj) 478 zvap = zvap * zcpy(ji,jj) 479 ENDIF 480 426 481 ! add to the general momentum trend 427 482 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap … … 445 500 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 446 501 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 502 503 IF(ln_wd) THEN 504 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 505 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 506 zuap = zuap * zcpx(ji,jj) 507 zvap = zvap * zcpy(ji,jj) 508 ENDIF 509 447 510 ! add to the general momentum trend 448 511 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap … … 467 530 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 468 531 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 532 533 IF(ln_wd) THEN 534 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 535 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 536 zuap = zuap * zcpx(ji,jj) 537 zvap = zvap * zcpy(ji,jj) 538 ENDIF 539 469 540 ! add to the general momentum trend 470 541 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap … … 476 547 ! 477 548 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 549 IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 478 550 ! 479 551 END SUBROUTINE hpg_sco … … 493 565 REAL(wp) :: z1_10, cffu, cffx ! " " 494 566 REAL(wp) :: z1_12, cffv, cffy ! " " 567 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 495 568 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 496 569 REAL(wp), POINTER, DIMENSION(:,:,:) :: dzx, dzy, dzz, dzu, dzv, dzw 497 570 REAL(wp), POINTER, DIMENSION(:,:,:) :: drhox, drhoy, drhoz, drhou, drhov, drhow 498 571 REAL(wp), POINTER, DIMENSION(:,:,:) :: rho_i, rho_j, rho_k 572 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter 499 573 !!---------------------------------------------------------------------- 500 574 ! … … 502 576 CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 503 577 CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 578 IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 504 579 ! 505 580 … … 508 583 IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend' 509 584 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, density Jacobian with cubic polynomial scheme' 585 ENDIF 586 587 IF(ln_wd) THEN 588 DO jj = 2, jpjm1 589 DO ji = 2, jpim1 590 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 591 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 592 & rn_wdmin1 + rn_wdmin2 593 594 IF(ll_tmp1) THEN 595 zcpx(ji,jj) = 1.0_wp 596 ELSE IF(ll_tmp2) THEN 597 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 598 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 599 & (sshn(ji+1,jj) - sshn(ji,jj))) 600 ELSE 601 zcpx(ji,jj) = 0._wp 602 END IF 603 604 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 605 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 606 & rn_wdmin1 + rn_wdmin2 607 608 IF(ll_tmp1) THEN 609 zcpy(ji,jj) = 1.0_wp 610 ELSE IF(ll_tmp2) THEN 611 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 612 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 613 & (sshn(ji,jj+1) - sshn(ji,jj))) 614 ELSE 615 zcpy(ji,jj) = 0._wp 616 END IF 617 END DO 618 END DO 619 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 510 620 ENDIF 511 621 … … 673 783 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) / e1u(ji,jj) 674 784 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) / e2v(ji,jj) 785 IF(ln_wd) THEN 786 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 787 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 788 ENDIF 675 789 ! add to the general momentum trend 676 790 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) … … 692 806 & + ( ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk ) ) & 693 807 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) / e2v(ji,jj) 808 IF(ln_wd) THEN 809 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 810 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 811 ENDIF 694 812 ! add to the general momentum trend 695 813 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) … … 702 820 CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 703 821 CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 822 IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 704 823 ! 705 824 END SUBROUTINE hpg_djc … … 727 846 !! The local variables for the correction term 728 847 INTEGER :: jk1, jis, jid, jjs, jjd 848 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 729 849 REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps 730 850 REAL(wp) :: zrhdt1 … … 732 852 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdept, zrhh 733 853 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 854 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter 734 855 !!---------------------------------------------------------------------- 735 856 ! 736 857 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 737 858 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 859 IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 738 860 ! 739 861 IF( kt == nit000 ) THEN … … 748 870 znad = 0.0_wp 749 871 IF( lk_vvl ) znad = 1._wp 872 873 IF(ln_wd) THEN 874 DO jj = 2, jpjm1 875 DO ji = 2, jpim1 876 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 877 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 878 & rn_wdmin1 + rn_wdmin2 879 880 IF(ll_tmp1) THEN 881 zcpx(ji,jj) = 1.0_wp 882 ELSE IF(ll_tmp2) THEN 883 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 884 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 885 & (sshn(ji+1,jj) - sshn(ji,jj))) 886 ELSE 887 zcpx(ji,jj) = 0._wp 888 END IF 889 890 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 891 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 892 & rn_wdmin1 + rn_wdmin2 893 894 IF(ll_tmp1.OR.ll_tmp2) THEN 895 zcpy(ji,jj) = 1.0_wp 896 ELSE IF(ll_tmp2) THEN 897 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 898 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 899 & (sshn(ji,jj+1) - sshn(ji,jj))) 900 ELSE 901 zcpy(ji,jj) = 0._wp 902 END IF 903 END DO 904 END DO 905 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 906 ENDIF 750 907 751 908 ! Clean 3-D work arrays … … 907 1064 ENDIF 908 1065 909 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * & 910 & umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 1066 IF(ln_wd) THEN 1067 zdpdx1 = zdpdx1 * zcpx(ji,jj) 1068 zdpdx2 = zdpdx2 * zcpx(ji,jj) 1069 ENDIF 1070 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) 911 1071 ENDIF 912 1072 … … 964 1124 ENDIF 965 1125 966 va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2)*& 967 & vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 1126 IF(ln_wd) THEN 1127 zdpdy1 = zdpdy1 * zcpy(ji,jj) 1128 zdpdy2 = zdpdy2 * zcpy(ji,jj) 1129 ENDIF 1130 1131 va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk) 968 1132 ENDIF 969 1133 … … 975 1139 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 976 1140 CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 1141 IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 977 1142 ! 978 1143 END SUBROUTINE hpg_prj -
branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r4624 r5066 11 11 !! 3.5 ! 2013-07 (J. Chanut) Switch to Forward-backward time stepping 12 12 !! 3.6 ! 2013-11 (A. Coward) Update for z-tilde compatibility 13 !! 3.? ! 2014-09 (H. Liu) Add Wetting/Drying components 13 14 !!--------------------------------------------------------------------- 14 15 #if defined key_dynspg_ts || defined key_esopa … … 40 41 USE timing ! Timing 41 42 USE sbcapr ! surface boundary condition: atmospheric pressure 43 USE wadlmt ! wetting/drying flux limter 42 44 USE dynadv, ONLY: ln_dynadv_vec 43 45 #if defined key_agrif … … 137 139 INTEGER, INTENT(in) :: kt ! ocean time-step index 138 140 ! 139 LOGICAL :: ll_fw_start ! if true, forward integration 140 LOGICAL :: ll_init ! if true, special startup of 2d equations 141 INTEGER :: ji, jj, jk, jn ! dummy loop indices 142 INTEGER :: ikbu, ikbv, noffset ! local integers 143 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 144 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 145 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - - 146 REAL(wp) :: zu_spg, zv_spg ! - - 147 REAL(wp) :: zhura, zhvra ! - - 148 REAL(wp) :: za0, za1, za2, za3 ! - - 141 LOGICAL :: ll_fw_start ! if true, forward integration 142 LOGICAL :: ll_init ! if true, special startup of 2d equations 143 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables used in W/D 144 INTEGER :: ji, jj, jk, jn ! dummy loop indices 145 INTEGER :: ikbu, ikbv, noffset ! local integers 146 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 147 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 148 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - - 149 REAL(wp) :: zu_spg, zv_spg ! - - 150 REAL(wp) :: zhura, zhvra ! - - 151 REAL(wp) :: za0, za1, za2, za3 ! - - 152 153 REAL(wp) :: ztmp ! temporary vaialbe 149 154 ! 150 155 REAL(wp), POINTER, DIMENSION(:,:) :: zun_e, zvn_e, zsshp2_e … … 154 159 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 155 160 REAL(wp), POINTER, DIMENSION(:,:) :: zhf 161 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef. 156 162 !!---------------------------------------------------------------------- 157 163 ! … … 165 171 CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a ) 166 172 CALL wrk_alloc( jpi, jpj, zhf ) 173 IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 167 174 ! 168 175 ! !* Local constant initialization … … 173 180 zraur = 1._wp / rau0 174 181 ! 175 IF( kt == nit000 .AND. neuler == 0 ) THEN 182 IF( kt == nit000 .AND. neuler == 0 ) THEN ! reciprocal of baroclinic time step 176 183 z2dt_bf = rdt 177 184 ELSE … … 180 187 z1_2dt_b = 1.0_wp / z2dt_bf 181 188 ! 182 ll_init = ln_bt_av 189 ll_init = ln_bt_av ! if no time averaging, then no specific restart 183 190 ll_fw_start = .FALSE. 184 191 ! 185 192 ! time offset in steps for bdy data update 186 193 IF (.NOT.ln_bt_fw) THEN ; noffset=-2*nn_baro ; ELSE ; noffset = 0 ; ENDIF 187 194 ! 188 IF( kt == nit000 ) THEN 195 IF( kt == nit000 ) THEN !* initialisation 189 196 ! 190 197 IF(lwp) WRITE(numout,*) … … 365 372 ! ! ---------------------------------------------------- 366 373 IF( lk_vvl ) THEN ! Variable volume : remove surface pressure gradient 367 DO jj = 2, jpjm1 368 DO ji = fs_2, fs_jpim1 ! vector opt. 369 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / e1u(ji,jj) 370 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / e2v(ji,jj) 371 END DO 372 END DO 374 IF(ln_wd) THEN ! Calculating and applying W/D gravity filters 375 DO jj = 2, jpjm1 376 DO ji = 2, jpim1 377 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 378 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 379 & rn_wdmin1 + rn_wdmin2 380 IF(ll_tmp1) THEN 381 zcpx(ji,jj) = 1.0_wp 382 ELSE IF(ll_tmp2) THEN 383 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 384 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 385 & (sshn(ji+1,jj) - sshn(ji,jj))) 386 ELSE 387 zcpx(ji,jj) = 0._wp 388 END IF 389 390 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 391 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 392 & rn_wdmin1 + rn_wdmin2 393 IF(ll_tmp1) THEN 394 zcpy(ji,jj) = 1.0_wp 395 ELSE IF(ll_tmp2) THEN 396 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 397 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) + bathy(ji,jj)) /& 398 & (sshn(ji,jj+1) - sshn(ji,jj))) 399 ELSE 400 zcpy(ji,jj) = 0._wp 401 END IF 402 END DO 403 END DO 404 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 405 ENDIF 406 407 IF(ln_wd) THEN ! Calculating and applying W/D gravity filters 408 DO jj = 2, jpjm1 409 DO ji = 2, jpim1 410 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / & 411 & e1u(ji,jj) * zcpx(ji,jj) 412 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / & 413 & e2v(ji,jj) * zcpy(ji,jj) 414 END DO 415 END DO 416 ELSE 417 DO jj = 2, jpjm1 418 DO ji = fs_2, fs_jpim1 ! vector opt. 419 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / e1u(ji,jj) 420 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / e2v(ji,jj) 421 END DO 422 END DO 423 END IF 373 424 ENDIF 374 425 … … 464 515 ! ! ==================== ! 465 516 ! Initialize barotropic variables: 517 518 IF(ln_wd) THEN !preserve the positivity of water depth 519 !ssh[b,n,a] should have already been processed for this 520 sshbb_e(:,:) = MAX(sshbb_e(:,:), rn_wdmin1 - bathy(:,:)) 521 sshb_e(:,:) = MAX(sshb_e(:,:) , rn_wdmin1 - bathy(:,:)) 522 ENDIF 523 466 524 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 467 525 sshn_e(:,:) = sshn (:,:) … … 537 595 zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points 538 596 zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 597 IF(ln_wd) THEN 598 zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1) 599 zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1) 600 END IF 539 601 ELSE 540 602 zhup2_e (:,:) = hu(:,:) … … 575 637 ENDIF 576 638 #endif 639 IF(ln_wd) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 577 640 ! 578 641 ! Sum over sub-time-steps to compute advective velocities … … 589 652 END DO 590 653 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1) 654 IF(ln_wd) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:)) 591 655 CALL lbc_lnk( ssha_e, 'T', 1._wp ) 592 656 … … 652 716 END DO 653 717 END DO 718 719 IF(ln_wd) THEN 720 zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 ) 721 zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 ) 722 END IF 723 !shall we call lbc_lnk for zhu(v)st_e() here? 724 654 725 ENDIF 655 726 ! … … 717 788 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:) 718 789 ! 719 ! Surface pressure trend: 720 DO jj = 2, jpjm1 721 DO ji = fs_2, fs_jpim1 ! vector opt. 722 ! Add surface pressure gradient 723 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 724 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 725 zwx(ji,jj) = zu_spg 726 zwy(ji,jj) = zv_spg 727 END DO 728 END DO 790 IF(ln_wd) THEN ! Calculating and applying W/D gravity filters 791 DO jj = 2, jpjm1 792 DO ji = 2, jpim1 793 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 794 ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 795 & rn_wdmin1 + rn_wdmin2 796 IF(ll_tmp1) THEN 797 zcpx(ji,jj) = 1.0_wp 798 ELSE IF(ll_tmp2) THEN 799 ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen ! here 800 zcpx(ji,jj) = ABS((zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) /& 801 & (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj))) 802 ELSE 803 zcpx(ji,jj) = 0._wp 804 END IF 805 806 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 807 ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 808 & rn_wdmin1 + rn_wdmin2 809 IF(ll_tmp1) THEN 810 zcpy(ji,jj) = 1.0_wp 811 ELSE IF(ll_tmp2) THEN 812 ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen ! here 813 zcpy(ji,jj) = ABS((zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) /& 814 & (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj))) 815 ELSE 816 zcpy(ji,jj) = 0._wp 817 END IF 818 END DO 819 END DO 820 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 821 ENDIF 822 823 IF(ln_wd) THEN 824 DO jj = 2, jpjm1 825 DO ji = 2, jpim1 826 ! Add surface pressure gradient 827 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 828 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 829 zwx(ji,jj) = zu_spg * zcpx(ji,jj) 830 zwy(ji,jj) = zv_spg * zcpy(ji,jj) 831 END DO 832 END DO 833 ELSE 834 DO jj = 2, jpjm1 835 DO ji = fs_2, fs_jpim1 ! vector opt. 836 ! Add surface pressure gradient 837 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 838 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 839 zwx(ji,jj) = zu_spg 840 zwy(ji,jj) = zv_spg 841 END DO 842 END DO 843 END IF 729 844 ! 730 845 ! Set next velocities: … … 750 865 DO ji = fs_2, fs_jpim1 ! vector opt. 751 866 752 zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1)) 753 zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 754 755 ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) & 756 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 867 IF(ln_wd) THEN 868 zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1) 869 zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1) 870 ELSE 871 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 872 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 873 END IF 874 875 zhura = umask(ji,jj,1)/(zhura + 1._wp - umask(ji,jj,1)) 876 zhvra = vmask(ji,jj,1)/(zhvra + 1._wp - vmask(ji,jj,1)) 877 878 ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) & 879 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 757 880 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & 758 881 & + hu(ji,jj) * zu_frc(ji,jj) ) & … … 770 893 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) 771 894 ! ! ---------------------------------------------- 772 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 773 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 895 IF(ln_wd) THEN 896 hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1) 897 hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1) 898 ELSE 899 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 900 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 901 END IF 902 774 903 hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 775 904 hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) … … 903 1032 CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a ) 904 1033 CALL wrk_dealloc( jpi, jpj, zhf ) 1034 IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 905 1035 ! 906 1036 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts') -
branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r4486 r5066 9 9 !! - ! 2010-09 (D.Storkey and E.O'Dea) bug fixes for BDY module 10 10 !! 3.3 ! 2011-10 (M. Leclair) split former ssh_wzv routine and remove all vvl related work 11 !! 3.? ! 2014-09 (H. Liu) add wetting and drying 11 12 !!---------------------------------------------------------------------- 12 13 … … 42 43 USE wrk_nemo ! Memory Allocation 43 44 USE timing ! Timing 45 USE wadlmt ! Wetting/Drying flux limting 44 46 45 47 IMPLICIT NONE … … 94 96 ENDIF 95 97 ! 96 CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity97 !98 98 z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog) 99 99 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 100 100 ! 101 z1_rau0 = 0.5_wp * r1_rau0 102 ! 103 IF(ln_wd) CALL wad_lmt(sshb, z1_rau0 * (emp_b(:,:) + emp(:,:)), z2dt) 104 105 CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity 106 ! 101 107 ! !------------------------------! 102 108 ! ! After Sea Surface Height ! … … 110 116 ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 111 117 ! 112 z1_rau0 = 0.5_wp * r1_rau0113 118 ssha(:,:) = ( sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * tmask(:,:,1) 114 119 -
branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
r4794 r5066 30 30 !! 3.4 ! 2011-11 (C. Harris) decomposition changes for running with CICE 31 31 !! ! 2012-05 (C. Calone, J. Simeon, G. Madec, C. Ethe) Add grid coarsening 32 !! 3.6.?! 2014-09 (H. Liu) Add Wetting and Drying 32 33 !!---------------------------------------------------------------------- 33 34 … … 233 234 NAMELIST/namcfg/ cp_cfg, cp_cfz, jp_cfg, jpidta, jpjdta, jpkdta, jpiglo, jpjglo, & 234 235 & jpizoom, jpjzoom, jperio 236 NAMELIST/namwad/ ln_wd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit 235 237 !!---------------------------------------------------------------------- 236 238 ! … … 257 259 READ ( numnam_cfg, namcfg, IOSTAT = ios, ERR = 904 ) 258 260 904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcfg in configuration namelist', .TRUE. ) 261 262 REWIND( numnam_ref ) ! Namelist namwad in reference namelist : Parameters for Wetting/Drying 263 READ ( numnam_ref, namwad, IOSTAT = ios, ERR = 905) 264 905 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in reference namelist', lwp ) 265 266 REWIND( numnam_cfg ) ! Namelist namwad in configuration namelist : Parameters for Wetting/Drying 267 READ ( numnam_cfg, namwad, IOSTAT = ios, ERR = 906) 268 906 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in configuration namelist', lwp ) 259 269 260 270 ! Force values for AGRIF zoom (cf. agrif_user.F90) … … 311 321 WRITE( numond, namctl ) 312 322 WRITE( numond, namcfg ) 323 WRITE( numond, namwad ) 313 324 ENDIF 314 325 … … 527 538 WRITE(numout,*) ' lateral cond. type (between 0 and 6) jperio = ', jperio 528 539 ENDIF 540 541 IF(lwp) THEN ! control print 542 WRITE(numout,*) 543 WRITE(numout,*) 'namwad : wetting and drying initialisation through namelist read' 544 WRITE(numout,*) '~~~~~~~ ' 545 WRITE(numout,*) ' Namelist namwad' 546 WRITE(numout,*) ' key to turn on/off wetting/drying ln_wd = ', ln_wd 547 WRITE(numout,*) ' minimum water depth on dried cells rn_wdmin1 = ', rn_wdmin1 548 WRITE(numout,*) ' tolerrance of min water depth on dried cells rn_wdmin2 = ', rn_wdmin2 549 WRITE(numout,*) ' land elevation below which wad considered rn_wdld = ', rn_wdld 550 WRITE(numout,*) ' maximum number of iteration for W/D limiter nn_wdit = ', nn_wdit 551 ENDIF 552 529 553 ! ! Parameter control 530 554 ! -
branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
r4328 r5066 118 118 USE trcstp ! passive tracer time-stepping (trc_stp routine) 119 119 #endif 120 USE wadlmt ! wetting/drying limiter (wad_lmt routine) 120 121 !!---------------------------------------------------------------------- 121 122 !! NEMO/OPA 3.3 , NEMO Consortium (2010)
Note: See TracChangeset
for help on using the changeset viewer.
