- Timestamp:
- 2012-07-11T13:22:58+02:00 (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r3226 r3432 54 54 ! ! ( rn_bb=0; top only, rn_bb =1; top and bottom) 55 55 REAL(wp) :: rn_hc = 150._wp ! Critical depth for s-sigma coordinates 56 56 PUBLIC rn_sbot_min, rn_sbot_max, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb,rn_hc 57 57 !! * Control permutation of array indices 58 58 # include "oce_ftrans.h90" … … 82 82 !! ln_zco=T z-coordinate 83 83 !! ln_zps=T z-coordinate with partial steps 84 !! ln_ sco=T s-coordinate84 !! ln_zco=T s-coordinate 85 85 !! 86 86 !! ** Action : define gdep., e3., mbathy and bathy … … 395 395 mbathy(:,:) = 0 ! set to zero extra halo points 396 396 bathy (:,:) = 0._wp ! (require for mpp case) 397 #if defined key_mpp_rkpart 398 DO jj = nldj, nlcj ! interior values 399 DO ji = nldi, nlci 400 #else 397 401 DO jj = 1, nlcj ! interior values 398 402 DO ji = 1, nlci 403 #endif 399 404 mbathy(ji,jj) = idta( mig(ji), mjg(jj) ) 400 405 bathy (ji,jj) = zdta( mig(ji), mjg(jj) ) … … 593 598 !! - update bathy : meter bathymetry (in meters) 594 599 !!---------------------------------------------------------------------- 595 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 596 USE wrk_nemo, ONLY: zbathy => wrk_2d_1 600 USE mapcomm_mod, ONLY: trimmed, nidx,eidx,sidx,widx 601 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 602 USE wrk_nemo, ONLY: zbathy => wrk_2d_1 597 603 !! 598 604 INTEGER :: ji, jj, jl ! dummy loop indices … … 636 642 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed' 637 643 ENDIF 644 638 645 IF( lk_mpp ) THEN 639 646 zbathy(:,:) = FLOAT( mbathy(:,:) ) … … 646 653 IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio 647 654 IF( lk_mpp ) THEN 648 IF( nbondi == -1 .OR. nbondi == 2) THEN655 IF( (nbondi == -1 .OR. nbondi == 2) .AND. (.NOT. trimmed(widx,narea) ) ) THEN 649 656 IF( jperio /= 1 ) mbathy(1,:) = 0 650 657 ENDIF 651 IF( nbondi == 1 .OR. nbondi == 2) THEN658 IF( (nbondi == 1 .OR. nbondi == 2) .AND. (.NOT. trimmed(eidx,narea) ) ) THEN 652 659 IF( jperio /= 1 ) mbathy(nlci,:) = 0 653 660 ENDIF … … 1188 1195 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 1189 1196 USE wrk_nemo, ONLY: zenv => wrk_2d_1 , ztmp => wrk_2d_2 , zmsk => wrk_2d_3 1190 USE wrk_nemo, ONLY: z ri => wrk_2d_4 , zrj => wrk_2d_5 , zhbat => wrk_2d_61197 USE wrk_nemo, ONLY: ztmp2 => wrk_2d_4 , zhbat => wrk_2d_5 1191 1198 USE wrk_nemo, ONLY: gsigw3 => wrk_3d_1 1192 1199 USE wrk_nemo, ONLY: gsigt3 => wrk_3d_2 … … 1199 1206 USE wrk_nemo, ONLY: esigwu3 => wrk_3d_9 1200 1207 USE wrk_nemo, ONLY: esigwv3 => wrk_3d_10 1208 USE mapcomm_mod, ONLY: trimmed, cyclic_bc 1209 USE mapcomm_mod, ONLY: nidx, eidx, sidx, widx 1210 ! USE arpdebugging, ONLY: dump_array 1201 1211 !! DCSE_NEMO: wrk_nemo module variables renamed, need additional directives 1202 1212 !FTRANS gsigw3 :I :I :z … … 1213 1223 INTEGER :: ji, jj, jk, jl ! dummy loop argument 1214 1224 INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers 1215 REAL(wp) :: zcoeft, zcoefw, zrmax, ztaper ! temporary scalars 1225 REAL(wp) :: zcoeft, zcoefw, zrmax, ztaper, zri, zrj ! temporary scalars 1226 REAL(wp), PARAMETER :: TOL_ZERO = 1.0E-20_wp ! Any value less than this assumed zero 1216 1227 ! 1217 1228 … … 1219 1230 !!---------------------------------------------------------------------- 1220 1231 1221 IF( wrk_in_use(2, 1,2,3,4,5 ,6) .OR. wrk_in_use(3, 1,2,3,4,5,6,7,8,9,10) ) THEN1232 IF( wrk_in_use(2, 1,2,3,4,5) .OR. wrk_in_use(3, 1,2,3,4,5,6,7,8,9,10) ) THEN 1222 1233 CALL ctl_stop('zgr_sco: ERROR - requested workspace arrays unavailable') ; RETURN 1223 1234 ENDIF … … 1269 1280 END DO 1270 1281 END DO 1282 1283 CALL lbc_lnk( zenv, 'T', 1._wp, lzero=.FALSE. ) 1271 1284 ! 1272 1285 ! Smooth the bathymetry (if required) … … 1282 1295 zrmax = 0._wp 1283 1296 zmsk(:,:) = 0._wp 1297 1284 1298 DO jj = 1, nlcj 1285 1299 DO ji = 1, nlci 1286 1300 iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi) 1287 ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj) 1288 zri(ji,jj) = ABS( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) 1289 zrj(ji,jj) = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) 1290 zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) ) 1291 IF( zri(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1._wp 1292 IF( zri(ji,jj) > rn_rmax ) zmsk(iip1,jj ) = 1._wp 1293 IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1._wp 1294 IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,ijp1) = 1._wp 1295 END DO 1296 END DO 1301 ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last row (jj=nclj+1 to jpj) 1302 zri = ABS( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) 1303 zrj = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) 1304 zrmax = MAX( zrmax, zri, zrj ) 1305 IF( zri > rn_rmax ) zmsk(ji ,jj ) = 1._wp 1306 IF( zri > rn_rmax ) zmsk(iip1,jj ) = 1._wp 1307 IF( zrj > rn_rmax ) zmsk(ji ,jj ) = 1._wp 1308 IF( zrj > rn_rmax ) zmsk(ji ,ijp1) = 1._wp 1309 END DO 1310 END DO 1311 1312 ! lateral boundary condition on zmsk: retain any 1's along closed 1313 ! boundary (use of lzero flag to lbc_lnk) 1314 CALL lbc_lnk( zmsk, 'T', 1._wp, lzero=.FALSE. ) 1315 1297 1316 IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain 1298 ! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX)1299 ztmp(:,:) = zmsk(:,:) ; CALL lbc_lnk( zmsk, 'T', 1._wp )1300 DO jj = 1, nlcj1301 DO ji = 1, nlci1302 zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )1303 END DO1304 END DO1305 1317 ! 1306 IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) ) 1318 IF(lwp)WRITE(numout,"('zgr_sco : iter=',I5,' rmax=',F8.4,' nb of pt= ',I8)") & 1319 jl, zrmax, INT( SUM(zmsk(:,:) ) ) 1307 1320 ! 1308 DO jj = 1, nlcj 1309 DO ji = 1, nlci 1321 !!$ IF(jl < 6)THEN ! .OR. (MOD(jl,1000) == 0) )THEN 1322 !!$ CALL dump_array(jl, 'zenv_before', zenv, withHalos=.TRUE.) 1323 !!$ CALL dump_array(jl, 'ztmp_before', ztmp, withHalos=.TRUE.) 1324 !!$ CALL dump_array(jl, 'zmsk_before', zmsk, withHalos=.TRUE.) 1325 !!$ END IF 1326 1327 ! Copy current surface before next smoothing iteration 1328 ztmp(:,:) = zenv(:,:) 1329 1330 DO jj = nldj, nlcj 1331 DO ji = nldi, nlci 1310 1332 iip1 = MIN( ji+1, nlci ) ! last line (ji=nlci) 1311 1333 ijp1 = MIN( jj+1, nlcj ) ! last raw (jj=nlcj) … … 1326 1348 END DO 1327 1349 ! 1328 DO jj = 1, nlcj 1329 DO ji = 1, nlci 1330 IF( zmsk(ji,jj) == 1._wp ) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) ) 1350 1351 ! Need to update halos of ztmp here but do not zero halos on closed 1352 ! boundaries 1353 CALL lbc_lnk( ztmp, 'T', 1._wp, lzero=.FALSE.) 1354 1355 DO jj = 1,nlcj 1356 DO ji = 1,nlci 1357 IF( zmsk(ji,jj) >= 1._wp-TOL_ZERO ) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) ) 1331 1358 END DO 1332 1359 END DO 1333 1360 ! 1334 ! Apply lateral boundary condition CAUTION: kept the value when the lbc field is zero 1335 ztmp(:,:) = zenv(:,:) ; CALL lbc_lnk( zenv, 'T', 1._wp ) 1336 DO jj = 1, nlcj 1337 DO ji = 1, nlci 1338 IF( zenv(ji,jj) == 0._wp ) zenv(ji,jj) = ztmp(ji,jj) 1339 END DO 1340 END DO 1361 ! Apply lateral boundary condition but do not zero on closed boundaries 1362 CALL lbc_lnk( zenv, 'T', 1._wp, lzero=.FALSE. ) 1363 1364 !!$ IF(jl < 6)THEN ! .OR. (MOD(jl,1000) == 0) )THEN 1365 !!$ CALL dump_array(jl, 'zenv', zenv, withHalos=.TRUE.) 1366 !!$ CALL dump_array(jl, 'ztmp', ztmp, withHalos=.TRUE.) 1367 !!$ CALL dump_array(jl, 'zmsk', zmsk, withHalos=.TRUE.) 1368 !!$ END IF 1369 1341 1370 ! ! ================ ! 1342 1371 END DO ! End loop ! … … 1365 1394 ENDIF 1366 1395 ENDIF 1396 1397 ! CALL dump_array(0, 'hbatt', hbatt, withHalos=.FALSE.) 1367 1398 1368 1399 ! ! ============================== … … 1576 1607 ! 1577 1608 !! H. Liu, POL. April 2009. Added for passing the scale check for the new released vvl code. 1609 where (e3t (:,:,:).eq.0.0) e3t(:,:,:) = 1.0 1610 where (e3u (:,:,:).eq.0.0) e3u(:,:,:) = 1.0 1611 where (e3v (:,:,:).eq.0.0) e3v(:,:,:) = 1.0 1612 where (e3f (:,:,:).eq.0.0) e3f(:,:,:) = 1.0 1613 where (e3w (:,:,:).eq.0.0) e3w(:,:,:) = 1.0 1614 where (e3uw (:,:,:).eq.0.0) e3uw(:,:,:) = 1.0 1615 where (e3vw (:,:,:).eq.0.0) e3vw(:,:,:) = 1.0 1616 1578 1617 1579 1618 fsdept(:,:,:) = gdept (:,:,:) … … 1597 1636 END DO 1598 1637 END DO 1638 1599 1639 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), & 1600 1640 & ' MAX ', MAXVAL( mbathy(:,:) ) … … 1642 1682 END DO 1643 1683 END DO 1644 DO jj = mj0(74), mj1(74) 1645 DO ji = mi0(100), mi1(100) 1646 WRITE(numout,*) 1647 WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj) 1648 WRITE(numout,*) ' ~~~~~~ --------------------' 1649 WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')") 1650 WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), & 1651 & fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk ) 1652 END DO 1653 END DO 1684 !!$ ARPDBG - out of bounds if jpj < 74 or jpi < 100, e.g. default GYRE 1685 !!$ DO jj = mj0(74), mj1(74) 1686 !!$ DO ji = mi0(100), mi1(100) 1687 !!$ WRITE(numout,*) 1688 !!$ WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj) 1689 !!$ WRITE(numout,*) ' ~~~~~~ --------------------' 1690 !!$ WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')") 1691 !!$ WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), & 1692 !!$ & fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk ) 1693 !!$ END DO 1694 !!$ END DO 1654 1695 ENDIF 1655 1696 … … 1668 1709 CALL ctl_stop( ctmp1 ) 1669 1710 ENDIF 1670 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN1711 IF( gdepw_1(ji,jj,jk) < 0._wp .OR. gdept_1(ji,jj,jk) < 0._wp ) THEN 1671 1712 WRITE(ctmp1,*) 'zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 1672 1713 CALL ctl_stop( ctmp1 ) … … 1677 1718 !!gm bug #endif 1678 1719 ! 1679 IF( wrk_not_released(2, 1,2,3,4,5 ,6) .OR. wrk_not_released(3, 1,2,3,4,5,6,7,8,9,10) ) &1720 IF( wrk_not_released(2, 1,2,3,4,5) .OR. wrk_not_released(3, 1,2,3,4,5,6,7,8,9,10) ) & 1680 1721 & CALL ctl_stop('dom:zgr_sco: failed to release workspace arrays') 1681 1722 !
Note: See TracChangeset
for help on using the changeset viewer.
