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 3926 – NEMO

Changeset 3926


Ignore:
Timestamp:
2013-06-17T14:50:15+02:00 (11 years ago)
Author:
jamesharle
Message:

Solution to ticket #1119: reproducibility using enveloping bathymetry.

File:
1 edited

Legend:

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

    r3764 r3926  
    169169      !!---------------------------------------------------------------------- 
    170170      !!                   ***  ROUTINE zgr_z  *** 
    171       !!                    
     171      !!                     
    172172      !! ** Purpose :   set the depth of model levels and the resulting  
    173173      !!      vertical scale factors. 
     
    639639         END DO 
    640640      END DO 
     641      IF( lk_mpp )   CALL mpp_sum( icompt ) 
    641642      IF( icompt == 0 ) THEN 
    642643         IF(lwp) WRITE(numout,*)'     no isolated ocean grid points' 
     
    11011102      INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers 
    11021103      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 
    11051109 
    11061110      NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & 
     
    11101114      IF( nn_timing == 1 )  CALL timing_start('zgr_sco') 
    11111115      ! 
    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     ! 
    11141119      REWIND( numnam )                       ! Read Namelist namzgr_sco : sigma-stretching parameters 
    11151120      READ  ( numnam, namzgr_sco ) 
     
    11581163      !                                        ! ============================= 
    11591164      ! 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(:,:) 
    11651172      !  
    11661173      ! Smooth the bathymetry (if required) 
     
    11701177      jl = 0 
    11711178      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         !                                                         ! ================ ! 
    11751195         jl = jl + 1 
    11761196         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 
    11781208         DO jj = 1, nlcj 
    11791209            DO ji = 1, nlci 
    11801210               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi) 
    11811211               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 
    11891222            END DO 
    11901223         END DO 
    11911224         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, nlcj 
    1195             DO ji = 1, nlci 
    1196                 zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) ) 
    1197             END DO 
    1198          END DO 
    11991225         ! 
    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 
    12011227         ! 
    12021228         DO jj = 1, nlcj 
    12031229            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) ) 
    12191231            END DO 
    12201232         END DO 
    12211233         ! 
    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 ) 
    12351235         !                                                  ! ================ ! 
    12361236      END DO                                                !     End loop     ! 
    12371237      !                                                     ! ================ ! 
    12381238      ! 
    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 
    12471244      ! 
    12481245      ! Envelope bathymetry saved in hbatt 
    12491246      hbatt(:,:) = zenv(:,:)  
     1247 
    12501248      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN 
    12511249         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' ) 
     
    14931491      END DO 
    14941492      ! 
    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                           )      ! 
    14971494      IF( nn_timing == 1 )  CALL timing_stop('zgr_sco') 
    14981495      ! 
Note: See TracChangeset for help on using the changeset viewer.