Changeset 15291

Timestamp:

2021-09-27T14:38:34+02:00 (3 years ago)

Author:

aumont

Message:

bug fix #2726

File:

• NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/TRP/trcsink.F90 (modified) (3 diffs)

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/TRP/trcsink.F90

@@ -79 +79 @@
                IF( tmask(ji,jj,jk) == 1.0 ) THEN
                    zwsmax =  0.5 * e3t(ji,jj,jk,Kmm) * rday / rsfact
-                   iiter(ji,jj) =  MAX( iiter(ji,jj), INT( pwsink(ji,jj,jk) / zwsmax ) )
+                   iiter(ji,jj) =  MAX( iiter(ji,jj), INT( pwsink(ji,jj,jk) / zwsmax ) + 1 )
                ENDIF
             END DO
@@ -126 +126 @@
       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   psinkflx  ! sinking fluxe
       !
-      INTEGER  ::   ji, jj, jk, jn
+      INTEGER  ::   ji, jj, jk, jn, jt
       REAL(wp) ::   zigma,zew,zign, zflx, zstep
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztraz, zakz, zwsink2, ztrb 
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztraz, zakz, zwsink2, ztrb, psinking 
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('trc_sink2')
       !
-      ztraz(:,:,:) = 0.e0
-      zakz (:,:,:) = 0.e0
-      ztrb (:,:,:) = tr(:,:,:,jp_tra,Kbb)
-
       DO jk = 1, jpkm1
          zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rday * tmask(:,:,jk+1) 
@@ -142 +138 @@
       zwsink2(:,:,1) = 0.e0
 
-
-      ! Vertical advective flux
-      DO jn = 1, 2
-         !  first guess of the slopes interior values
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-            !
-            zstep = rsfact / REAL( kiter(ji,jj), wp ) / 2.
-            !              
-            DO jk = 2, jpkm1
-               ztraz(ji,jj,jk) = ( tr(ji,jj,jk-1,jp_tra,Kbb) - tr(ji,jj,jk,jp_tra,Kbb) ) * tmask(ji,jj,jk)
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         ! Vertical advective flux
+         zstep = rsfact / REAL( kiter(ji,jj), wp ) / 2.
+         DO jt = 1, kiter(ji,jj)
+            ztraz(ji,jj,:) = 0.e0
+            zakz (ji,jj,:) = 0.e0
+            ztrb (ji,jj,:) = tr(ji,jj,:,jp_tra,Kbb)
+            DO jn = 1, 2
+               !              
+               DO jk = 2, jpkm1
+                  ztraz(ji,jj,jk) = ( tr(ji,jj,jk-1,jp_tra,Kbb) - tr(ji,jj,jk,jp_tra,Kbb) ) * tmask(ji,jj,jk)
+               END DO
+               ztraz(ji,jj,1  ) = 0.0
+               ztraz(ji,jj,jpk) = 0.0
+
+               ! slopes
+               DO jk = 2, jpkm1
+                  zign = 0.25 + SIGN( 0.25_wp, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) )
+                  zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign
+               END DO
+      
+               ! Slopes limitation
+               DO jk = 2, jpkm1
+                  zakz(ji,jj,jk) = SIGN( 1.0_wp, zakz(ji,jj,jk) ) *        &
+                     &             MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) )
+               END DO
+      
+               ! vertical advective flux
+               DO jk = 1, jpkm1
+                  zigma = zwsink2(ji,jj,jk+1) * zstep / e3w(ji,jj,jk+1,Kmm)
+                  zew   = zwsink2(ji,jj,jk+1)
+                  psinking(ji,jj,jk+1) = -zew * ( tr(ji,jj,jk,jp_tra,Kbb) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep
+               END DO
+               !
+               ! Boundary conditions
+               psinking(ji,jj,1  ) = 0.e0
+               psinking(ji,jj,jpk) = 0.e0
+      
+               DO jk = 1, jpkm1
+                  zflx = ( psinking(ji,jj,jk) - psinking(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm)
+                  tr(ji,jj,jk,jp_tra,Kbb) = tr(ji,jj,jk,jp_tra,Kbb) + zflx
+               END DO
             END DO
-            ztraz(ji,jj,1  ) = 0.0
-            ztraz(ji,jj,jpk) = 0.0
-
-            ! slopes
-            DO jk = 2, jpkm1
-               zign = 0.25 + SIGN( 0.25_wp, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) )
-               zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign
+            DO jk = 1, jpkm1
+               zflx = ( psinking(ji,jj,jk) - psinking(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm)
+               ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx
             END DO
-      
-            ! Slopes limitation
-            DO jk = 2, jpkm1
-               zakz(ji,jj,jk) = SIGN( 1.0_wp, zakz(ji,jj,jk) ) *        &
-                  &             MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) )
-            END DO
-      
-            ! vertical advective flux
-            DO jk = 1, jpkm1
-               zigma = zwsink2(ji,jj,jk+1) * zstep / e3w(ji,jj,jk+1,Kmm)
-               zew   = zwsink2(ji,jj,jk+1)
-               psinkflx(ji,jj,jk+1) = -zew * ( tr(ji,jj,jk,jp_tra,Kbb) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep
-            END DO
-            !
-            ! Boundary conditions
-            psinkflx(ji,jj,1  ) = 0.e0
-            psinkflx(ji,jj,jpk) = 0.e0
-      
-            DO jk=1,jpkm1
-               zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm)
-               tr(ji,jj,jk,jp_tra,Kbb) = tr(ji,jj,jk,jp_tra,Kbb) + zflx
-            END DO
-         END_2D
-      END DO
-
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
-         zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm)
-         ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx
-      END_3D
-
-      tr(:,:,:,jp_tra,Kbb) = ztrb(:,:,:)
-      psinkflx(:,:,:)   = 2. * psinkflx(:,:,:)
+
+            tr(ji,jj,:,jp_tra,Kbb) = ztrb(ji,jj,:)
+            psinkflx(ji,jj,:)   = psinkflx(ji,jj,:) + 2. * psinking(ji,jj,:)
+         END DO
+      END_2D
       !
       IF( ln_timing )  CALL timing_stop('trc_sink2')