- Timestamp:
- 2013-11-05T13:25:45+01:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r4148 r4153 1106 1106 INTEGER :: ios ! Local integer output status for namelist read 1107 1107 REAL(wp) :: zrmax, ztaper ! temporary scalars 1108 REAL(wp) :: zrfact ! temporary scalars 1109 REAL(wp), POINTER, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 1110 1111 ! 1112 REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, zri, zrj, zhbat 1108 ! 1109 REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 1113 1110 1114 1111 NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & … … 1118 1115 IF( nn_timing == 1 ) CALL timing_start('zgr_sco') 1119 1116 ! 1120 CALL wrk_alloc( jpi, jpj, ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) 1121 CALL wrk_alloc( jpi, jpj, zenv, zri, zrj, zhbat ) 1122 ! 1117 CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat ) 1118 ! 1123 1119 REWIND( numnam_ref ) ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters 1124 1120 READ ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901) … … 1173 1169 ! ! ============================= 1174 1170 ! use r-value to create hybrid coordinates 1175 ! DO jj = 1, jpj 1176 ! DO ji = 1, jpi 1177 ! zenv(ji,jj) = MAX( bathy(ji,jj), 0._wp ) 1178 ! END DO 1179 ! END DO 1180 ! CALL lbc_lnk( zenv, 'T', 1._wp ) 1181 zenv(:,:) = bathy(:,:) 1171 DO jj = 1, jpj 1172 DO ji = 1, jpi 1173 zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min ) 1174 END DO 1175 END DO 1182 1176 ! 1183 1177 ! Smooth the bathymetry (if required) … … 1187 1181 jl = 0 1188 1182 zrmax = 1._wp 1189 ! 1190 ! set scaling factor used in reducing vertical gradients 1191 zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax ) 1192 ! 1193 ! initialise temporary evelope depth arrays 1194 ztmpi1(:,:) = zenv(:,:) 1195 ztmpi2(:,:) = zenv(:,:) 1196 ztmpj1(:,:) = zenv(:,:) 1197 ztmpj2(:,:) = zenv(:,:) 1198 ! 1199 ! initialise temporary r-value arrays 1200 zri(:,:) = 1._wp 1201 zrj(:,:) = 1._wp 1202 ! ! ================ ! 1203 DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) ! Iterative loop ! 1204 ! ! ================ ! 1183 ! ! ================ ! 1184 DO WHILE( jl <= 10000 .AND. zrmax > rn_rmax ) ! Iterative loop ! 1185 ! ! ================ ! 1205 1186 jl = jl + 1 1206 1187 zrmax = 0._wp 1207 ! we set zrmax from previous r-values (zri abd zrj) first 1208 ! if set after current r-value calculation (as previously) 1209 ! we could exit DO WHILE prematurely before checking r-value 1210 ! of current zenv 1211 DO jj = 1, nlcj 1212 DO ji = 1, nlci 1213 zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) ) 1214 END DO 1215 END DO 1216 zri(:,:) = 0._wp 1217 zrj(:,:) = 0._wp 1188 zmsk(:,:) = 0._wp 1218 1189 DO jj = 1, nlcj 1219 1190 DO ji = 1, nlci 1220 1191 iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi) 1221 1192 ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj) 1222 IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN 1223 zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) 1224 END IF 1225 IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN 1226 zrj(ji,jj) = ( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) 1227 END IF 1228 IF( zri(ji,jj) > rn_rmax ) ztmpi1(ji ,jj ) = zenv(iip1,jj ) * zrfact 1229 IF( zri(ji,jj) < -rn_rmax ) ztmpi2(iip1,jj ) = zenv(ji ,jj ) * zrfact 1230 IF( zrj(ji,jj) > rn_rmax ) ztmpj1(ji ,jj ) = zenv(ji ,ijp1) * zrfact 1231 IF( zrj(ji,jj) < -rn_rmax ) ztmpj2(ji ,ijp1) = zenv(ji ,jj ) * zrfact 1193 zri(ji,jj) = ABS( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) 1194 zrj(ji,jj) = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) 1195 zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) ) 1196 IF( zri(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1._wp 1197 IF( zri(ji,jj) > rn_rmax ) zmsk(iip1,jj ) = 1._wp 1198 IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1._wp 1199 IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,ijp1) = 1._wp 1232 1200 END DO 1233 1201 END DO 1234 1202 IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain 1203 ! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX) 1204 ztmp(:,:) = zmsk(:,:) ; CALL lbc_lnk( zmsk, 'T', 1._wp ) 1205 DO jj = 1, nlcj 1206 DO ji = 1, nlci 1207 zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) ) 1208 END DO 1209 END DO 1235 1210 ! 1236 IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax 1211 IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) ) 1237 1212 ! 1238 1213 DO jj = 1, nlcj 1239 1214 DO ji = 1, nlci 1240 zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) ) 1215 iip1 = MIN( ji+1, nlci ) ! last line (ji=nlci) 1216 ijp1 = MIN( jj+1, nlcj ) ! last raw (jj=nlcj) 1217 iim1 = MAX( ji-1, 1 ) ! first line (ji=nlci) 1218 ijm1 = MAX( jj-1, 1 ) ! first raw (jj=nlcj) 1219 IF( zmsk(ji,jj) == 1._wp ) THEN 1220 ztmp(ji,jj) = ( & 1221 & zenv(iim1,ijp1)*zmsk(iim1,ijp1) + zenv(ji,ijp1)*zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1) & 1222 & + zenv(iim1,jj )*zmsk(iim1,jj ) + zenv(ji,jj )* 2._wp + zenv(iip1,jj )*zmsk(iip1,jj ) & 1223 & + zenv(iim1,ijm1)*zmsk(iim1,ijm1) + zenv(ji,ijm1)*zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1) & 1224 & ) / ( & 1225 & zmsk(iim1,ijp1) + zmsk(ji,ijp1) + zmsk(iip1,ijp1) & 1226 & + zmsk(iim1,jj ) + 2._wp + zmsk(iip1,jj ) & 1227 & + zmsk(iim1,ijm1) + zmsk(ji,ijm1) + zmsk(iip1,ijm1) & 1228 & ) 1229 ENDIF 1241 1230 END DO 1242 1231 END DO 1243 1232 ! 1244 CALL lbc_lnk( zenv, 'T', 1._wp ) 1233 DO jj = 1, nlcj 1234 DO ji = 1, nlci 1235 IF( zmsk(ji,jj) == 1._wp ) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) ) 1236 END DO 1237 END DO 1238 ! 1239 ! Apply lateral boundary condition CAUTION: keep the value when the lbc field is zero 1240 ztmp(:,:) = zenv(:,:) ; CALL lbc_lnk( zenv, 'T', 1._wp ) 1241 DO jj = 1, nlcj 1242 DO ji = 1, nlci 1243 IF( zenv(ji,jj) == 0._wp ) zenv(ji,jj) = ztmp(ji,jj) 1244 END DO 1245 END DO 1245 1246 ! ! ================ ! 1246 1247 END DO ! End loop ! 1247 1248 ! ! ================ ! 1248 1249 ! 1249 ! DO jj = 1, jpj 1250 ! DO ji = 1, jpi 1251 ! zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale values 1252 ! END DO 1253 ! END DO 1250 ! Fill ghost rows with appropriate values to avoid undefined e3 values with some mpp decompositions 1251 DO ji = nlci+1, jpi 1252 zenv(ji,1:nlcj) = zenv(nlci,1:nlcj) 1253 END DO 1254 ! 1255 DO jj = nlcj+1, jpj 1256 zenv(:,jj) = zenv(:,nlcj) 1257 END DO 1254 1258 ! 1255 1259 ! Envelope bathymetry saved in hbatt 1256 1260 hbatt(:,:) = zenv(:,:) 1257 1258 1261 IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN 1259 1262 CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' ) 1260 1263 DO jj = 1, jpj 1261 1264 DO ji = 1, jpi 1262 ztaper = EXP( -(gphit(ji,jj)/8._wp)**2 )1265 ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp ) 1263 1266 hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper ) 1264 1267 END DO … … 1375 1378 fsde3w(:,:,:) = gdep3w(:,:,:) 1376 1379 ! 1377 where (e3t (:,:,:).eq.0.0) e3t(:,:,:) = 1.0 1378 where (e3u (:,:,:).eq.0.0) e3u(:,:,:) = 1.0 1379 where (e3v (:,:,:).eq.0.0) e3v(:,:,:) = 1.0 1380 where (e3f (:,:,:).eq.0.0) e3f(:,:,:) = 1.0 1381 where (e3w (:,:,:).eq.0.0) e3w(:,:,:) = 1.0 1382 where (e3uw (:,:,:).eq.0.0) e3uw(:,:,:) = 1.0 1383 where (e3vw (:,:,:).eq.0.0) e3vw(:,:,:) = 1.0 1384 1380 where (e3t (:,:,:).eq.0.0) e3t(:,:,:) = 1._wp 1381 where (e3u (:,:,:).eq.0.0) e3u(:,:,:) = 1._wp 1382 where (e3v (:,:,:).eq.0.0) e3v(:,:,:) = 1._wp 1383 where (e3f (:,:,:).eq.0.0) e3f(:,:,:) = 1._wp 1384 where (e3w (:,:,:).eq.0.0) e3w(:,:,:) = 1._wp 1385 where (e3uw (:,:,:).eq.0.0) e3uw(:,:,:) = 1._wp 1386 where (e3vw (:,:,:).eq.0.0) e3vw(:,:,:) = 1._wp 1387 1388 #if defined key_agrif 1389 ! Ensure meaningful vertical scale factors in ghost lines/columns 1390 IF( .NOT. Agrif_Root() ) THEN 1391 ! 1392 IF((nbondi == -1).OR.(nbondi == 2)) THEN 1393 e3u(1,:,:) = e3u(2,:,:) 1394 ENDIF 1395 ! 1396 IF((nbondi == 1).OR.(nbondi == 2)) THEN 1397 e3u(nlci-1,:,:) = e3u(nlci-2,:,:) 1398 ENDIF 1399 ! 1400 IF((nbondj == -1).OR.(nbondj == 2)) THEN 1401 e3v(:,1,:) = e3v(:,2,:) 1402 ENDIF 1403 ! 1404 IF((nbondj == 1).OR.(nbondj == 2)) THEN 1405 e3v(:,nlcj-1,:) = e3v(:,nlcj-2,:) 1406 ENDIF 1407 ! 1408 ENDIF 1409 #endif 1385 1410 1386 1411 fsdept(:,:,:) = gdept (:,:,:) … … 1431 1456 WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(1,1,jk), fsdepw(1,1,jk), & 1432 1457 & fse3t (1,1,jk), fse3w (1,1,jk), jk=1,jpk ) 1433 DO jj = mj0(20), mj1(20) 1434 DO ji = mi0(20), mi1(20) 1458 iip1 = MIN(20, jpiglo-1) ! for config with i smaller than 20 points 1459 ijp1 = MIN(20, jpjglo-1) ! for config with j smaller than 20 points 1460 DO jj = mj0(ijp1), mj1(ijp1) 1461 DO ji = mi0(iip1), mi1(iip1) 1435 1462 WRITE(numout,*) 1436 WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj) 1463 WRITE(numout,*) ' domzgr: vertical coordinates : point (',iip1,',',ijp1,',k) bathy = ', & 1464 & bathy(ji,jj), hbatt(ji,jj) 1437 1465 WRITE(numout,*) ' ~~~~~~ --------------------' 1438 1466 WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')") … … 1441 1469 END DO 1442 1470 END DO 1443 DO jj = mj0(74), mj1(74) 1444 DO ji = mi0(100), mi1(100) 1471 iip1 = MIN( 74, jpiglo-1) 1472 ijp1 = MIN( 100, jpjglo-1) 1473 DO jj = mj0(ijp1), mj1(ijp1) 1474 DO ji = mi0(iip1), mi1(iip1) 1445 1475 WRITE(numout,*) 1446 WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj) 1476 WRITE(numout,*) ' domzgr: vertical coordinates : point (',iip1,',',ijp1,',k) bathy = ', & 1477 & bathy(ji,jj), hbatt(ji,jj) 1447 1478 WRITE(numout,*) ' ~~~~~~ --------------------' 1448 1479 WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')") … … 1501 1532 END DO 1502 1533 ! 1503 CALL wrk_dealloc( jpi, jpj, zenv, ztmpi1, ztmpi2, ztmpj1, ztmpj2, zri, zrj, zhbat ) ! 1534 CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat ) 1535 ! 1504 1536 IF( nn_timing == 1 ) CALL timing_stop('zgr_sco') 1505 1537 ! … … 1730 1762 ENDDO 1731 1763 ! 1732 CALL lbc_lnk(e3t ,'T',1.) ; CALL lbc_lnk(e3u ,'T',1.)1733 CALL lbc_lnk(e3v ,'T',1.) ; CALL lbc_lnk(e3f ,'T',1.)1734 CALL lbc_lnk(e3w ,'T',1.)1735 CALL lbc_lnk(e3uw,'T',1.) ; CALL lbc_lnk(e3vw,'T',1.)1736 !1737 1764 ! ! ============= 1738 1765 … … 1831 1858 !!---------------------------------------------------------------------- 1832 1859 ! 1833 pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1 ) + rn_thetb ) ) &1860 pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1,wp) + rn_thetb ) ) & 1834 1861 & - TANH( rn_thetb * rn_theta ) ) & 1835 1862 & * ( COSH( rn_theta ) & … … 1857 1884 ! 1858 1885 IF ( rn_theta == 0 ) then ! uniform sigma 1859 pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )1886 pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1,wp ) 1860 1887 ELSE ! stretched sigma 1861 pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1 )) ) ) / SINH( rn_theta ) &1862 & + pbb * ( (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1 )) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta ) ) &1888 pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1,wp)) ) ) / SINH( rn_theta ) & 1889 & + pbb * ( (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1,wp)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta ) ) & 1863 1890 & / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) ) 1864 1891 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
