- Timestamp:
- 2013-11-04T13:54:28+01:00 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r4147 r4148 176 176 !!---------------------------------------------------------------------- 177 177 !! *** ROUTINE zgr_z *** 178 !! 178 !! 179 179 !! ** Purpose : set the depth of model levels and the resulting 180 180 !! vertical scale factors. … … 645 645 END DO 646 646 END DO 647 IF( lk_mpp ) CALL mpp_sum( icompt ) 647 648 IF( icompt == 0 ) THEN 648 649 IF(lwp) WRITE(numout,*)' no isolated ocean grid points' … … 1105 1106 INTEGER :: ios ! Local integer output status for namelist read 1106 1107 REAL(wp) :: zrmax, ztaper ! temporary scalars 1107 ! 1108 REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 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 1109 1113 1110 1114 NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & … … 1114 1118 IF( nn_timing == 1 ) CALL timing_start('zgr_sco') 1115 1119 ! 1116 CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat ) 1117 ! 1120 CALL wrk_alloc( jpi, jpj, ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) 1121 CALL wrk_alloc( jpi, jpj, zenv, zri, zrj, zhbat ) 1122 ! 1118 1123 REWIND( numnam_ref ) ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters 1119 1124 READ ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901) … … 1168 1173 ! ! ============================= 1169 1174 ! use r-value to create hybrid coordinates 1170 DO jj = 1, jpj 1171 DO ji = 1, jpi 1172 zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min ) 1173 END DO 1174 END DO 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(:,:) 1175 1182 ! 1176 1183 ! Smooth the bathymetry (if required) … … 1180 1187 jl = 0 1181 1188 zrmax = 1._wp 1182 ! ! ================ ! 1183 DO WHILE( jl <= 10000 .AND. zrmax > rn_rmax ) ! Iterative loop ! 1184 ! ! ================ ! 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 ! ! ================ ! 1185 1205 jl = jl + 1 1186 1206 zrmax = 0._wp 1187 zmsk(:,:) = 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 1218 DO jj = 1, nlcj 1189 1219 DO ji = 1, nlci 1190 1220 iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi) 1191 1221 ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj) 1192 zri(ji,jj) = ABS( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) 1193 zrj(ji,jj) = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) 1194 zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) ) 1195 IF( zri(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1._wp 1196 IF( zri(ji,jj) > rn_rmax ) zmsk(iip1,jj ) = 1._wp 1197 IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1._wp 1198 IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,ijp1) = 1._wp 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 1199 1232 END DO 1200 1233 END DO 1201 1234 IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain 1202 ! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX)1203 ztmp(:,:) = zmsk(:,:) ; CALL lbc_lnk( zmsk, 'T', 1._wp )1204 DO jj = 1, nlcj1205 DO ji = 1, nlci1206 zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )1207 END DO1208 END DO1209 1235 ! 1210 IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax , ' nb of pt= ', INT( SUM(zmsk(:,:) ) )1236 IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax 1211 1237 ! 1212 1238 DO jj = 1, nlcj 1213 1239 DO ji = 1, nlci 1214 iip1 = MIN( ji+1, nlci ) ! last line (ji=nlci) 1215 ijp1 = MIN( jj+1, nlcj ) ! last raw (jj=nlcj) 1216 iim1 = MAX( ji-1, 1 ) ! first line (ji=nlci) 1217 ijm1 = MAX( jj-1, 1 ) ! first raw (jj=nlcj) 1218 IF( zmsk(ji,jj) == 1._wp ) THEN 1219 ztmp(ji,jj) = ( & 1220 & zenv(iim1,ijp1)*zmsk(iim1,ijp1) + zenv(ji,ijp1)*zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1) & 1221 & + zenv(iim1,jj )*zmsk(iim1,jj ) + zenv(ji,jj )* 2._wp + zenv(iip1,jj )*zmsk(iip1,jj ) & 1222 & + zenv(iim1,ijm1)*zmsk(iim1,ijm1) + zenv(ji,ijm1)*zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1) & 1223 & ) / ( & 1224 & zmsk(iim1,ijp1) + zmsk(ji,ijp1) + zmsk(iip1,ijp1) & 1225 & + zmsk(iim1,jj ) + 2._wp + zmsk(iip1,jj ) & 1226 & + zmsk(iim1,ijm1) + zmsk(ji,ijm1) + zmsk(iip1,ijm1) & 1227 & ) 1228 ENDIF 1240 zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) ) 1229 1241 END DO 1230 1242 END DO 1231 1243 ! 1232 DO jj = 1, nlcj 1233 DO ji = 1, nlci 1234 IF( zmsk(ji,jj) == 1._wp ) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) ) 1235 END DO 1236 END DO 1237 ! 1238 ! Apply lateral boundary condition CAUTION: keep the value when the lbc field is zero 1239 ztmp(:,:) = zenv(:,:) ; CALL lbc_lnk( zenv, 'T', 1._wp ) 1240 DO jj = 1, nlcj 1241 DO ji = 1, nlci 1242 IF( zenv(ji,jj) == 0._wp ) zenv(ji,jj) = ztmp(ji,jj) 1243 END DO 1244 END DO 1244 CALL lbc_lnk( zenv, 'T', 1._wp ) 1245 1245 ! ! ================ ! 1246 1246 END DO ! End loop ! 1247 1247 ! ! ================ ! 1248 1248 ! 1249 ! Fill ghost rows with appropriate values to avoid undefined e3 values with some mpp decompositions 1250 DO ji = nlci+1, jpi 1251 zenv(ji,1:nlcj) = zenv(nlci,1:nlcj) 1252 END DO 1253 ! 1254 DO jj = nlcj+1, jpj 1255 zenv(:,jj) = zenv(:,nlcj) 1256 END DO 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 1257 1254 ! 1258 1255 ! Envelope bathymetry saved in hbatt 1259 1256 hbatt(:,:) = zenv(:,:) 1257 1260 1258 IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN 1261 1259 CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' ) … … 1503 1501 END DO 1504 1502 ! 1505 CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat ) 1506 ! 1503 CALL wrk_dealloc( jpi, jpj, zenv, ztmpi1, ztmpi2, ztmpj1, ztmpj2, zri, zrj, zhbat ) ! 1507 1504 IF( nn_timing == 1 ) CALL timing_stop('zgr_sco') 1508 1505 !
Note: See TracChangeset
for help on using the changeset viewer.
