Changeset 4245 for branches/2013/dev_LOCEAN_CMCC_INGV_MERC_UKMO_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
- Timestamp:
- 2013-11-19T12:19:21+01:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_LOCEAN_CMCC_INGV_MERC_UKMO_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r4153 r4245 374 374 IF(lwp) WRITE(numout,*) 375 375 IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin' 376 idta(:,:) = jpkm1 ! before last level 377 zdta(:,:) = gdepw_0(jpk) ! last w-point depth 378 h_oce = gdepw_0(jpk) 376 IF( rn_bathy > 0.01 ) THEN 377 IF(lwp) WRITE(numout,*) ' Depth = rn_bathy read in namelist' 378 zdta(:,:) = rn_bathy 379 IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk 380 idta(:,:) = jpkm1 381 ELSE ! z-coordinate (zco or zps): step-like topography 382 idta(:,:) = jpkm1 383 DO jk = 1, jpkm1 384 WHERE( gdept_0(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_0(jk+1) ) idta(:,:) = jk 385 END DO 386 ENDIF 387 ELSE 388 IF(lwp) WRITE(numout,*) ' Depth = depthw(jpkm1)' 389 idta(:,:) = jpkm1 ! before last level 390 zdta(:,:) = gdepw_0(jpk) ! last w-point depth 391 h_oce = gdepw_0(jpk) 392 ENDIF 379 393 ELSE ! bump centered in the basin 380 394 IF(lwp) WRITE(numout,*) … … 1106 1120 INTEGER :: ios ! Local integer output status for namelist read 1107 1121 REAL(wp) :: zrmax, ztaper ! temporary scalars 1108 ! 1122 REAL(wp) :: zrfact 1123 ! 1124 REAL(wp), POINTER, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 1109 1125 REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 1110 1126 … … 1115 1131 IF( nn_timing == 1 ) CALL timing_start('zgr_sco') 1116 1132 ! 1117 CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat)1133 CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) 1118 1134 ! 1119 1135 REWIND( numnam_ref ) ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters … … 1169 1185 ! ! ============================= 1170 1186 ! use r-value to create hybrid coordinates 1187 zenv(:,:) = bathy(:,:) 1188 ! 1189 ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing 1171 1190 DO jj = 1, jpj 1172 1191 DO ji = 1, jpi 1173 zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min ) 1174 END DO 1175 END DO 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 1200 ENDIF 1201 ENDIF 1202 END DO 1203 END DO 1204 ! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero 1205 CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' ) 1176 1206 ! 1177 ! Smooth the bathymetry (if required)1207 ! smooth the bathymetry (if required) 1178 1208 scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea) 1179 1209 scobot(:,:) = bathy(:,:) ! ocean bottom depth … … 1181 1211 jl = 0 1182 1212 zrmax = 1._wp 1183 ! ! ================ ! 1184 DO WHILE( jl <= 10000 .AND. zrmax > rn_rmax ) ! Iterative loop ! 1185 ! ! ================ ! 1213 ! 1214 ! 1215 ! set scaling factor used in reducing vertical gradients 1216 zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax ) 1217 ! 1218 ! initialise temporary evelope depth arrays 1219 ztmpi1(:,:) = zenv(:,:) 1220 ztmpi2(:,:) = zenv(:,:) 1221 ztmpj1(:,:) = zenv(:,:) 1222 ztmpj2(:,:) = zenv(:,:) 1223 ! 1224 ! initialise temporary r-value arrays 1225 zri(:,:) = 1._wp 1226 zrj(:,:) = 1._wp 1227 ! ! ================ ! 1228 DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) ! Iterative loop ! 1229 ! ! ================ ! 1186 1230 jl = jl + 1 1187 1231 zrmax = 0._wp 1188 zmsk(:,:) = 0._wp 1232 ! we set zrmax from previous r-values (zri and zrj) first 1233 ! if set after current r-value calculation (as previously) 1234 ! we could exit DO WHILE prematurely before checking r-value 1235 ! of current zenv 1236 DO jj = 1, nlcj 1237 DO ji = 1, nlci 1238 zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) ) 1239 END DO 1240 END DO 1241 zri(:,:) = 0._wp 1242 zrj(:,:) = 0._wp 1189 1243 DO jj = 1, nlcj 1190 1244 DO ji = 1, nlci 1191 1245 iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi) 1192 1246 ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj) 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 1247 IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN 1248 zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) 1249 END IF 1250 IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN 1251 zrj(ji,jj) = ( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) 1252 END IF 1253 IF( zri(ji,jj) > rn_rmax ) ztmpi1(ji ,jj ) = zenv(iip1,jj ) * zrfact 1254 IF( zri(ji,jj) < -rn_rmax ) ztmpi2(iip1,jj ) = zenv(ji ,jj ) * zrfact 1255 IF( zrj(ji,jj) > rn_rmax ) ztmpj1(ji ,jj ) = zenv(ji ,ijp1) * zrfact 1256 IF( zrj(ji,jj) < -rn_rmax ) ztmpj2(ji ,ijp1) = zenv(ji ,jj ) * zrfact 1200 1257 END DO 1201 1258 END DO 1202 1259 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, nlcj1206 DO ji = 1, nlci1207 zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )1208 END DO1209 END DO1210 1260 ! 1211 IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax , ' nb of pt= ', INT( SUM(zmsk(:,:) ) )1261 IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax 1212 1262 ! 1213 1263 DO jj = 1, nlcj 1214 1264 DO ji = 1, nlci 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 1230 END DO 1231 END DO 1232 ! 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 1265 zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) ) 1266 END DO 1267 END DO 1268 ! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero 1269 CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' ) 1246 1270 ! ! ================ ! 1247 1271 END DO ! End loop ! 1248 1272 ! ! ================ ! 1249 ! 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) 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 1257 1277 END DO 1258 1278 ! … … 1532 1552 END DO 1533 1553 ! 1534 CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat)1554 CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) 1535 1555 ! 1536 1556 IF( nn_timing == 1 ) CALL timing_stop('zgr_sco')
Note: See TracChangeset
for help on using the changeset viewer.
