Changeset 1073 for trunk/NEMO/TOP_SRC/PISCES/p4zsink.F90

Timestamp:

2008-06-05T14:15:34+02:00 (16 years ago)

Author:

cetlod

Message:

update PISCES model, see ticket:190

File:

• trunk/NEMO/TOP_SRC/PISCES/p4zsink.F90 (modified) (2 diffs)

trunk/NEMO/TOP_SRC/PISCES/p4zsink.F90

@@ -13 +13 @@
    USE oce_trc         !
    USE trp_trc
-   USE sms
-   USE p4zsink2        !
+   USE sms_pisces
    USE prtctl_trc
 
@@ -615 +614 @@
 #endif
 
+   SUBROUTINE p4z_sink2( pwsink, psinkflx, jp_tra )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_sink2  ***
+      !!
+      !! ** Purpose :   Compute the sedimentation terms for the various sinking
+      !!     particles. The scheme used to compute the trends is based
+      !!     on MUSCL.
+      !!
+      !! ** Method  : - this ROUTINE compute not exactly the advection but the
+      !!      transport term, i.e.  div(u*tra).
+      !!---------------------------------------------------------------------
+      INTEGER , INTENT(in   )                         ::   jp_tra    ! tracer index index      
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwsink    ! sinking speed
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   psinkflx  ! sinking fluxe
+      !!
+      INTEGER  ::   ji, jj, jk, jn
+      REAL(wp) ::   zigma,zew,zstep,zign, zflx
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  ztraz, zakz
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zwsink2
+      !!---------------------------------------------------------------------
+
+      zstep  = rfact2 / 2.
+
+      ztraz(:,:,:) = 0.e0
+      zakz (:,:,:) = 0.e0
+
+      DO jk = 1, jpkm1
+# if defined key_off_degrad
+         zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rjjss * tmask(:,:,jk+1) * facvol(:,:,jk)
+# else
+         zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rjjss * tmask(:,:,jk+1)
+# endif
+      END DO
+      zwsink2(:,:,1) = 0.e0
+
+
+      ! Vertical advective flux
+      DO jn = 1, 2
+         !  first guess of the slopes interior values
+         DO jk = 2, jpkm1
+            ztraz(:,:,jk) = ( trn(:,:,jk-1,jp_tra) - trn(:,:,jk,jp_tra) ) * tmask(:,:,jk)
+         END DO
+         ztraz(:,:,1  ) = 0.0
+         ztraz(:,:,jpk) = 0.0
+
+         ! slopes
+         DO jk = 2, jpkm1
+            DO jj = 1,jpj
+               DO ji = 1, jpi
+                  zign = 0.25 + SIGN( 0.25, 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
+            END DO
+         END DO
+         
+         ! Slopes limitation
+         DO jk = 2, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zakz(ji,jj,jk) = SIGN( 1., 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
+            END DO
+         END DO
+         
+         ! vertical advective flux
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj      
+               DO ji = 1, jpi    
+                  zigma = zwsink2(ji,jj,jk+1) * zstep / fse3w(ji,jj,jk+1)
+                  zew   = zwsink2(ji,jj,jk+1)
+                  psinkflx(ji,jj,jk+1) = -zew * ( trn(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep
+               END DO
+            END DO
+         END DO
+         !
+         ! Boundary conditions
+         psinkflx(:,:,1  ) = 0.e0
+         psinkflx(:,:,jpk) = 0.e0
+         
+         DO jk=1,jpkm1
+            DO jj = 1,jpj
+               DO ji = 1, jpi
+                  zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
+                  trn(ji,jj,jk,jp_tra) = trn(ji,jj,jk,jp_tra) + zflx
+               END DO
+            END DO
+         END DO
+
+      ENDDO
+
+      DO jk=1,jpkm1
+         DO jj = 1,jpj
+            DO ji = 1, jpi
+               zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
+               trb(ji,jj,jk,jp_tra) = trb(ji,jj,jk,jp_tra) + 2. * zflx
+            END DO
+         END DO
+      END DO
+
+      trn(:,:,:,jp_tra) = trb(:,:,:,jp_tra)
+      psinkflx(:,:,:)   = 2. * psinkflx(:,:,:)
+
+      !
+   END SUBROUTINE p4z_sink2
+
 #else
    !!======================================================================