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 15325 for NEMO/trunk/src/TOP/TRP/trcsink.F90 – NEMO

Ignore:
Timestamp:
2021-10-04T14:07:47+02:00 (3 years ago)
Author:
cetlod
Message:

bugfix in TOP/TRP/trcsink.F90, see ticket #2726

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/TOP/TRP/trcsink.F90

    r15090 r15325  
    7979               IF( tmask(ji,jj,jk) == 1.0 ) THEN 
    8080                   zwsmax =  0.5 * e3t(ji,jj,jk,Kmm) * rday / rsfact 
    81                    iiter(ji,jj) =  MAX( iiter(ji,jj), INT( pwsink(ji,jj,jk) / zwsmax ) ) 
     81                   iiter(ji,jj) =  MAX( iiter(ji,jj), INT( pwsink(ji,jj,jk) / zwsmax ) + 1 ) 
    8282               ENDIF 
    8383            END DO 
     
    126126      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   psinkflx  ! sinking fluxe 
    127127      ! 
    128       INTEGER  ::   ji, jj, jk, jn 
     128      INTEGER  ::   ji, jj, jk, jn, jt 
    129129      REAL(wp) ::   zigma,zew,zign, zflx, zstep 
    130       REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztraz, zakz, zwsink2, ztrb  
     130      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztraz, zakz, zwsink2, ztrb, psinking  
    131131      !!--------------------------------------------------------------------- 
    132132      ! 
    133133      IF( ln_timing )   CALL timing_start('trc_sink2') 
    134134      ! 
    135       ztraz(:,:,:) = 0.e0 
    136       zakz (:,:,:) = 0.e0 
    137       ztrb (:,:,:) = tr(:,:,:,jp_tra,Kbb) 
    138  
    139135      DO jk = 1, jpkm1 
    140136         zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rday * tmask(:,:,jk+1)  
     
    142138      zwsink2(:,:,1) = 0.e0 
    143139 
    144  
    145       ! Vertical advective flux 
    146       DO jn = 1, 2 
    147          !  first guess of the slopes interior values 
    148          DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    149             ! 
    150             zstep = rsfact / REAL( kiter(ji,jj), wp ) / 2. 
    151             !               
    152             DO jk = 2, jpkm1 
    153                ztraz(ji,jj,jk) = ( tr(ji,jj,jk-1,jp_tra,Kbb) - tr(ji,jj,jk,jp_tra,Kbb) ) * tmask(ji,jj,jk) 
     140      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     141         ! Vertical advective flux 
     142         zstep = rsfact / REAL( kiter(ji,jj), wp ) / 2. 
     143         DO jt = 1, kiter(ji,jj) 
     144            ztraz(ji,jj,:) = 0.e0 
     145            zakz (ji,jj,:) = 0.e0 
     146            ztrb (ji,jj,:) = tr(ji,jj,:,jp_tra,Kbb) 
     147            DO jn = 1, 2 
     148               !               
     149               DO jk = 2, jpkm1 
     150                  ztraz(ji,jj,jk) = ( tr(ji,jj,jk-1,jp_tra,Kbb) - tr(ji,jj,jk,jp_tra,Kbb) ) * tmask(ji,jj,jk) 
     151               END DO 
     152               ztraz(ji,jj,1  ) = 0.0 
     153               ztraz(ji,jj,jpk) = 0.0 
     154 
     155               ! slopes 
     156               DO jk = 2, jpkm1 
     157                  zign = 0.25 + SIGN( 0.25_wp, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) ) 
     158                  zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign 
     159               END DO 
     160       
     161               ! Slopes limitation 
     162               DO jk = 2, jpkm1 
     163                  zakz(ji,jj,jk) = SIGN( 1.0_wp, zakz(ji,jj,jk) ) *        & 
     164                     &             MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) ) 
     165               END DO 
     166       
     167               ! vertical advective flux 
     168               DO jk = 1, jpkm1 
     169                  zigma = zwsink2(ji,jj,jk+1) * zstep / e3w(ji,jj,jk+1,Kmm) 
     170                  zew   = zwsink2(ji,jj,jk+1) 
     171                  psinking(ji,jj,jk+1) = -zew * ( tr(ji,jj,jk,jp_tra,Kbb) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep 
     172               END DO 
     173               ! 
     174               ! Boundary conditions 
     175               psinking(ji,jj,1  ) = 0.e0 
     176               psinking(ji,jj,jpk) = 0.e0 
     177       
     178               DO jk = 1, jpkm1 
     179                  zflx = ( psinking(ji,jj,jk) - psinking(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm) 
     180                  tr(ji,jj,jk,jp_tra,Kbb) = tr(ji,jj,jk,jp_tra,Kbb) + zflx 
     181               END DO 
    154182            END DO 
    155             ztraz(ji,jj,1  ) = 0.0 
    156             ztraz(ji,jj,jpk) = 0.0 
    157  
    158             ! slopes 
    159             DO jk = 2, jpkm1 
    160                zign = 0.25 + SIGN( 0.25_wp, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) ) 
    161                zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign 
     183            DO jk = 1, jpkm1 
     184               zflx = ( psinking(ji,jj,jk) - psinking(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm) 
     185               ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx 
    162186            END DO 
    163        
    164             ! Slopes limitation 
    165             DO jk = 2, jpkm1 
    166                zakz(ji,jj,jk) = SIGN( 1.0_wp, zakz(ji,jj,jk) ) *        & 
    167                   &             MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) ) 
    168             END DO 
    169        
    170             ! vertical advective flux 
    171             DO jk = 1, jpkm1 
    172                zigma = zwsink2(ji,jj,jk+1) * zstep / e3w(ji,jj,jk+1,Kmm) 
    173                zew   = zwsink2(ji,jj,jk+1) 
    174                psinkflx(ji,jj,jk+1) = -zew * ( tr(ji,jj,jk,jp_tra,Kbb) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep 
    175             END DO 
    176             ! 
    177             ! Boundary conditions 
    178             psinkflx(ji,jj,1  ) = 0.e0 
    179             psinkflx(ji,jj,jpk) = 0.e0 
    180        
    181             DO jk=1,jpkm1 
    182                zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm) 
    183                tr(ji,jj,jk,jp_tra,Kbb) = tr(ji,jj,jk,jp_tra,Kbb) + zflx 
    184             END DO 
    185          END_2D 
    186       END DO 
    187  
    188       DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 
    189          zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm) 
    190          ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx 
    191       END_3D 
    192  
    193       tr(:,:,:,jp_tra,Kbb) = ztrb(:,:,:) 
    194       psinkflx(:,:,:)   = 2. * psinkflx(:,:,:) 
     187 
     188            tr(ji,jj,:,jp_tra,Kbb) = ztrb(ji,jj,:) 
     189            psinkflx(ji,jj,:)   = psinkflx(ji,jj,:) + 2. * psinking(ji,jj,:) 
     190         END DO 
     191      END_2D 
    195192      ! 
    196193      IF( ln_timing )  CALL timing_stop('trc_sink2') 
Note: See TracChangeset for help on using the changeset viewer.