Changeset 3998 for trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
- Timestamp:
- 2013-08-02T19:31:01+02:00 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r3971 r3998 1102 1102 INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers 1103 1103 REAL(wp) :: zrmax, ztaper ! temporary scalars 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 1104 ! 1105 REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 1109 1106 1110 1107 NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & … … 1114 1111 IF( nn_timing == 1 ) CALL timing_start('zgr_sco') 1115 1112 ! 1116 CALL wrk_alloc( jpi, jpj, ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) 1117 CALL wrk_alloc( jpi, jpj, zenv, zri, zrj, zhbat ) 1118 ! 1113 CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat ) 1114 ! 1119 1115 REWIND( numnam ) ! Read Namelist namzgr_sco : sigma-stretching parameters 1120 1116 READ ( numnam, namzgr_sco ) … … 1163 1159 ! ! ============================= 1164 1160 ! use r-value to create hybrid coordinates 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(:,:) 1161 DO jj = 1, jpj 1162 DO ji = 1, jpi 1163 zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min ) 1164 END DO 1165 END DO 1172 1166 ! 1173 1167 ! Smooth the bathymetry (if required) … … 1177 1171 jl = 0 1178 1172 zrmax = 1._wp 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 ! ! ================ ! 1173 ! ! ================ ! 1174 DO WHILE( jl <= 10000 .AND. zrmax > rn_rmax ) ! Iterative loop ! 1175 ! ! ================ ! 1195 1176 jl = jl + 1 1196 1177 zrmax = 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 zmsk(:,:) = 0._wp 1208 1179 DO jj = 1, nlcj 1209 1180 DO ji = 1, nlci 1210 1181 iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi) 1211 1182 ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj) 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 1183 zri(ji,jj) = ABS( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) 1184 zrj(ji,jj) = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) 1185 zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) ) 1186 IF( zri(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1._wp 1187 IF( zri(ji,jj) > rn_rmax ) zmsk(iip1,jj ) = 1._wp 1188 IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1._wp 1189 IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,ijp1) = 1._wp 1222 1190 END DO 1223 1191 END DO 1224 1192 IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain 1193 ! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX) 1194 ztmp(:,:) = zmsk(:,:) ; CALL lbc_lnk( zmsk, 'T', 1._wp ) 1195 DO jj = 1, nlcj 1196 DO ji = 1, nlci 1197 zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) ) 1198 END DO 1199 END DO 1225 1200 ! 1226 IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax 1201 IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) ) 1227 1202 ! 1228 1203 DO jj = 1, nlcj 1229 1204 DO ji = 1, nlci 1230 zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) ) 1205 iip1 = MIN( ji+1, nlci ) ! last line (ji=nlci) 1206 ijp1 = MIN( jj+1, nlcj ) ! last raw (jj=nlcj) 1207 iim1 = MAX( ji-1, 1 ) ! first line (ji=nlci) 1208 ijm1 = MAX( jj-1, 1 ) ! first raw (jj=nlcj) 1209 IF( zmsk(ji,jj) == 1._wp ) THEN 1210 ztmp(ji,jj) = ( & 1211 & zenv(iim1,ijp1)*zmsk(iim1,ijp1) + zenv(ji,ijp1)*zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1) & 1212 & + zenv(iim1,jj )*zmsk(iim1,jj ) + zenv(ji,jj )* 2._wp + zenv(iip1,jj )*zmsk(iip1,jj ) & 1213 & + zenv(iim1,ijm1)*zmsk(iim1,ijm1) + zenv(ji,ijm1)*zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1) & 1214 & ) / ( & 1215 & zmsk(iim1,ijp1) + zmsk(ji,ijp1) + zmsk(iip1,ijp1) & 1216 & + zmsk(iim1,jj ) + 2._wp + zmsk(iip1,jj ) & 1217 & + zmsk(iim1,ijm1) + zmsk(ji,ijm1) + zmsk(iip1,ijm1) & 1218 & ) 1219 ENDIF 1231 1220 END DO 1232 1221 END DO 1233 1222 ! 1234 CALL lbc_lnk( zenv, 'T', 1._wp ) 1223 DO jj = 1, nlcj 1224 DO ji = 1, nlci 1225 IF( zmsk(ji,jj) == 1._wp ) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) ) 1226 END DO 1227 END DO 1228 ! 1229 ! Apply lateral boundary condition CAUTION: keep the value when the lbc field is zero 1230 ztmp(:,:) = zenv(:,:) ; CALL lbc_lnk( zenv, 'T', 1._wp ) 1231 DO jj = 1, nlcj 1232 DO ji = 1, nlci 1233 IF( zenv(ji,jj) == 0._wp ) zenv(ji,jj) = ztmp(ji,jj) 1234 END DO 1235 END DO 1235 1236 ! ! ================ ! 1236 1237 END DO ! End loop ! 1237 1238 ! ! ================ ! 1238 1239 ! 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 1240 ! Fill ghost rows with appropriate values to avoid undefined e3 values with some mpp decompositions 1241 DO ji = nlci+1, jpi 1242 zenv(ji,1:nlcj) = zenv(nlci,1:nlcj) 1243 END DO 1244 ! 1245 DO jj = nlcj+1, jpj 1246 zenv(:,jj) = zenv(:,nlcj) 1247 END DO 1244 1248 ! 1245 1249 ! Envelope bathymetry saved in hbatt 1246 1250 hbatt(:,:) = zenv(:,:) 1247 1248 1251 IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN 1249 1252 CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' ) … … 1519 1522 END DO 1520 1523 ! 1521 CALL wrk_dealloc( jpi, jpj, zenv, ztmpi1, ztmpi2, ztmpj1, ztmpj2, zri, zrj, zhbat ) ! 1524 CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat ) 1525 ! 1522 1526 IF( nn_timing == 1 ) CALL timing_stop('zgr_sco') 1523 1527 !
Note: See TracChangeset
for help on using the changeset viewer.