Changeset 3926 for trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
- Timestamp:
- 2013-06-17T14:50:15+02:00 (10 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r3764 r3926 169 169 !!---------------------------------------------------------------------- 170 170 !! *** ROUTINE zgr_z *** 171 !! 171 !! 172 172 !! ** Purpose : set the depth of model levels and the resulting 173 173 !! vertical scale factors. … … 639 639 END DO 640 640 END DO 641 IF( lk_mpp ) CALL mpp_sum( icompt ) 641 642 IF( icompt == 0 ) THEN 642 643 IF(lwp) WRITE(numout,*)' no isolated ocean grid points' … … 1101 1102 INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers 1102 1103 REAL(wp) :: zrmax, ztaper ! temporary scalars 1103 ! 1104 REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 1104 REAL(wp) :: zrfact ! temporary scalars 1105 REAL(wp), POINTER, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 1106 1107 ! 1108 REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, zri, zrj, zhbat 1105 1109 1106 1110 NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & … … 1110 1114 IF( nn_timing == 1 ) CALL timing_start('zgr_sco') 1111 1115 ! 1112 CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat ) 1113 ! 1116 CALL wrk_alloc( jpi, jpj, ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) 1117 CALL wrk_alloc( jpi, jpj, zenv, zri, zrj, zhbat ) 1118 ! 1114 1119 REWIND( numnam ) ! Read Namelist namzgr_sco : sigma-stretching parameters 1115 1120 READ ( numnam, namzgr_sco ) … … 1158 1163 ! ! ============================= 1159 1164 ! use r-value to create hybrid coordinates 1160 DO jj = 1, jpj 1161 DO ji = 1, jpi 1162 zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min ) 1163 END DO 1164 END DO 1165 ! DO jj = 1, jpj 1166 ! DO ji = 1, jpi 1167 ! zenv(ji,jj) = MAX( bathy(ji,jj), 0._wp ) 1168 ! END DO 1169 ! END DO 1170 ! CALL lbc_lnk( zenv, 'T', 1._wp ) 1171 zenv(:,:) = bathy(:,:) 1165 1172 ! 1166 1173 ! Smooth the bathymetry (if required) … … 1170 1177 jl = 0 1171 1178 zrmax = 1._wp 1172 ! ! ================ ! 1173 DO WHILE( jl <= 10000 .AND. zrmax > rn_rmax ) ! Iterative loop ! 1174 ! ! ================ ! 1179 ! 1180 ! set scaling factor used in reducing vertical gradients 1181 zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax ) 1182 ! 1183 ! initialise temporary evelope depth arrays 1184 ztmpi1(:,:) = zenv(:,:) 1185 ztmpi2(:,:) = zenv(:,:) 1186 ztmpj1(:,:) = zenv(:,:) 1187 ztmpj2(:,:) = zenv(:,:) 1188 ! 1189 ! initialise temporary r-value arrays 1190 zri(:,:) = 1._wp 1191 zrj(:,:) = 1._wp 1192 ! ! ================ ! 1193 DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) ! Iterative loop ! 1194 ! ! ================ ! 1175 1195 jl = jl + 1 1176 1196 zrmax = 0._wp 1177 zmsk(:,:) = 0._wp 1197 ! we set zrmax from previous r-values (zri abd zrj) first 1198 ! if set after current r-value calculation (as previously) 1199 ! we could exit DO WHILE prematurely before checking r-value 1200 ! of current zenv 1201 DO jj = 1, nlcj 1202 DO ji = 1, nlci 1203 zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) ) 1204 END DO 1205 END DO 1206 zri(:,:) = 0._wp 1207 zrj(:,:) = 0._wp 1178 1208 DO jj = 1, nlcj 1179 1209 DO ji = 1, nlci 1180 1210 iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi) 1181 1211 ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj) 1182 zri(ji,jj) = ABS( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) 1183 zrj(ji,jj) = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) 1184 zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) ) 1185 IF( zri(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1._wp 1186 IF( zri(ji,jj) > rn_rmax ) zmsk(iip1,jj ) = 1._wp 1187 IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1._wp 1188 IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,ijp1) = 1._wp 1212 IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN 1213 zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) 1214 END IF 1215 IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN 1216 zrj(ji,jj) = ( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) 1217 END IF 1218 IF( zri(ji,jj) > rn_rmax ) ztmpi1(ji ,jj ) = zenv(iip1,jj ) * zrfact 1219 IF( zri(ji,jj) < -rn_rmax ) ztmpi2(iip1,jj ) = zenv(ji ,jj ) * zrfact 1220 IF( zrj(ji,jj) > rn_rmax ) ztmpj1(ji ,jj ) = zenv(ji ,ijp1) * zrfact 1221 IF( zrj(ji,jj) < -rn_rmax ) ztmpj2(ji ,ijp1) = zenv(ji ,jj ) * zrfact 1189 1222 END DO 1190 1223 END DO 1191 1224 IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain 1192 ! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX)1193 ztmp(:,:) = zmsk(:,:) ; CALL lbc_lnk( zmsk, 'T', 1._wp )1194 DO jj = 1, nlcj1195 DO ji = 1, nlci1196 zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )1197 END DO1198 END DO1199 1225 ! 1200 IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax , ' nb of pt= ', INT( SUM(zmsk(:,:) ) )1226 IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax 1201 1227 ! 1202 1228 DO jj = 1, nlcj 1203 1229 DO ji = 1, nlci 1204 iip1 = MIN( ji+1, nlci ) ! last line (ji=nlci) 1205 ijp1 = MIN( jj+1, nlcj ) ! last raw (jj=nlcj) 1206 iim1 = MAX( ji-1, 1 ) ! first line (ji=nlci) 1207 ijm1 = MAX( jj-1, 1 ) ! first raw (jj=nlcj) 1208 IF( zmsk(ji,jj) == 1._wp ) THEN 1209 ztmp(ji,jj) = ( & 1210 & zenv(iim1,ijp1)*zmsk(iim1,ijp1) + zenv(ji,ijp1)*zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1) & 1211 & + zenv(iim1,jj )*zmsk(iim1,jj ) + zenv(ji,jj )* 2._wp + zenv(iip1,jj )*zmsk(iip1,jj ) & 1212 & + zenv(iim1,ijm1)*zmsk(iim1,ijm1) + zenv(ji,ijm1)*zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1) & 1213 & ) / ( & 1214 & zmsk(iim1,ijp1) + zmsk(ji,ijp1) + zmsk(iip1,ijp1) & 1215 & + zmsk(iim1,jj ) + 2._wp + zmsk(iip1,jj ) & 1216 & + zmsk(iim1,ijm1) + zmsk(ji,ijm1) + zmsk(iip1,ijm1) & 1217 & ) 1218 ENDIF 1230 zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) ) 1219 1231 END DO 1220 1232 END DO 1221 1233 ! 1222 DO jj = 1, nlcj 1223 DO ji = 1, nlci 1224 IF( zmsk(ji,jj) == 1._wp ) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) ) 1225 END DO 1226 END DO 1227 ! 1228 ! Apply lateral boundary condition CAUTION: keep the value when the lbc field is zero 1229 ztmp(:,:) = zenv(:,:) ; CALL lbc_lnk( zenv, 'T', 1._wp ) 1230 DO jj = 1, nlcj 1231 DO ji = 1, nlci 1232 IF( zenv(ji,jj) == 0._wp ) zenv(ji,jj) = ztmp(ji,jj) 1233 END DO 1234 END DO 1234 CALL lbc_lnk( zenv, 'T', 1._wp ) 1235 1235 ! ! ================ ! 1236 1236 END DO ! End loop ! 1237 1237 ! ! ================ ! 1238 1238 ! 1239 ! Fill ghost rows with appropriate values to avoid undefined e3 values with some mpp decompositions 1240 DO ji = nlci+1, jpi 1241 zenv(ji,1:nlcj) = zenv(nlci,1:nlcj) 1242 END DO 1243 ! 1244 DO jj = nlcj+1, jpj 1245 zenv(:,jj) = zenv(:,nlcj) 1246 END DO 1239 ! DO jj = 1, jpj 1240 ! DO ji = 1, jpi 1241 ! zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale values 1242 ! END DO 1243 ! END DO 1247 1244 ! 1248 1245 ! Envelope bathymetry saved in hbatt 1249 1246 hbatt(:,:) = zenv(:,:) 1247 1250 1248 IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN 1251 1249 CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' ) … … 1493 1491 END DO 1494 1492 ! 1495 CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat ) 1496 ! 1493 CALL wrk_dealloc( jpi, jpj, zenv, ztmpi1, ztmpi2, ztmpj1, ztmpj2, zri, zrj, zhbat ) ! 1497 1494 IF( nn_timing == 1 ) CALL timing_stop('zgr_sco') 1498 1495 !
Note: See TracChangeset
for help on using the changeset viewer.
