New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 3998 for trunk – NEMO

Changeset 3998 for trunk


Ignore:
Timestamp:
2013-08-02T19:31:01+02:00 (11 years ago)
Author:
jamesharle
Message:

remove new evelope smoothing code and revert back to r3925 version of the code - see ticket 1119

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r3971 r3998  
    11021102      INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers 
    11031103      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 
    11091106 
    11101107      NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & 
     
    11141111      IF( nn_timing == 1 )  CALL timing_start('zgr_sco') 
    11151112      ! 
    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      ! 
    11191115      REWIND( numnam )                       ! Read Namelist namzgr_sco : sigma-stretching parameters 
    11201116      READ  ( numnam, namzgr_sco ) 
     
    11631159      !                                        ! ============================= 
    11641160      ! 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 
    11721166      !  
    11731167      ! Smooth the bathymetry (if required) 
     
    11771171      jl = 0 
    11781172      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         !                                                  ! ================ ! 
    11951176         jl = jl + 1 
    11961177         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 
    12081179         DO jj = 1, nlcj 
    12091180            DO ji = 1, nlci 
    12101181               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi) 
    12111182               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 
    12221190            END DO 
    12231191         END DO 
    12241192         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 
    12251200         ! 
    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(:,:) ) ) 
    12271202         ! 
    12281203         DO jj = 1, nlcj 
    12291204            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 
    12311220            END DO 
    12321221         END DO 
    12331222         ! 
    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 
    12351236         !                                                  ! ================ ! 
    12361237      END DO                                                !     End loop     ! 
    12371238      !                                                     ! ================ ! 
    12381239      ! 
    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 
    12441248      ! 
    12451249      ! Envelope bathymetry saved in hbatt 
    12461250      hbatt(:,:) = zenv(:,:)  
    1247  
    12481251      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN 
    12491252         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' ) 
     
    15191522      END DO 
    15201523      ! 
    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      ! 
    15221526      IF( nn_timing == 1 )  CALL timing_stop('zgr_sco') 
    15231527      ! 
Note: See TracChangeset for help on using the changeset viewer.