Changeset 2082

Timestamp:

2010-09-10T12:32:58+02:00 (14 years ago)

Author:

cetlod

Message:

Improve the merge of TRA-TRC, see ticket #717

Location:

branches/DEV_r2006_merge_TRA_TRC/NEMO

Files:

• NST_SRC/agrif_user.F90 (modified) (2 diffs)
• OFF_SRC/dtadyn.F90 (modified) (7 diffs)
• OFF_SRC/istate.F90 (modified) (1 diff)
• OFF_SRC/opa.F90 (modified) (2 diffs)
• OPA_SRC/DIA/diaar5.F90 (modified) (3 diffs)
• OPA_SRC/IOM/restart.F90 (modified) (2 diffs)
• OPA_SRC/TRA/eosbn2.F90 (modified) (22 diffs)
• OPA_SRC/TRA/traadv.F90 (modified) (4 diffs)
• OPA_SRC/TRA/traadv_cen2.F90 (modified) (4 diffs)
• OPA_SRC/TRA/traadv_eiv.F90 (modified) (2 diffs)
• OPA_SRC/TRA/traadv_muscl.F90 (modified) (7 diffs)
• OPA_SRC/TRA/traadv_muscl2.F90 (modified) (7 diffs)
• OPA_SRC/TRA/traadv_qck.F90 (modified) (10 diffs)
• OPA_SRC/TRA/traadv_tvd.F90 (modified) (11 diffs)
• OPA_SRC/TRA/traadv_ubs.F90 (modified) (11 diffs)
• OPA_SRC/TRA/trabbl.F90 (modified) (3 diffs)
• OPA_SRC/TRA/traldf.F90 (modified) (2 diffs)
• OPA_SRC/TRA/traldf_bilap.F90 (modified) (2 diffs)
• OPA_SRC/TRA/traldf_bilapg.F90 (modified) (2 diffs)
• OPA_SRC/TRA/traldf_iso.F90 (modified) (2 diffs)
• OPA_SRC/TRA/traldf_lap.F90 (modified) (2 diffs)
• OPA_SRC/TRA/tranpc.F90 (modified) (1 diff)
• OPA_SRC/TRA/tranxt.F90 (modified) (4 diffs)
• OPA_SRC/TRA/trazdf.F90 (modified) (1 diff)
• OPA_SRC/TRA/trazdf_exp.F90 (modified) (2 diffs)
• OPA_SRC/TRA/trazdf_imp.F90 (modified) (2 diffs)
• OPA_SRC/TRA/zpshde.F90 (modified) (7 diffs)
• OPA_SRC/TRD/trdicp.F90 (modified) (1 diff)
• OPA_SRC/istate.F90 (modified) (2 diffs)
• OPA_SRC/oce.F90 (modified) (1 diff)
• OPA_SRC/opa.F90 (modified) (8 diffs)
• OPA_SRC/step.F90 (modified) (4 diffs)
• OPA_SRC/step_oce.F90 (modified) (1 diff)
• OPA_SRC/trc_oce.F90 (modified) (1 diff)
• TOP_SRC/C14b/trcsms_c14b.F90 (modified) (2 diffs)
• TOP_SRC/CFC/trcsms_cfc.F90 (modified) (2 diffs)
• TOP_SRC/LOBSTER/trcbio.F90 (modified) (1 diff)
• TOP_SRC/LOBSTER/trcexp.F90 (modified) (1 diff)
• TOP_SRC/LOBSTER/trcopt.F90 (modified) (1 diff)
• TOP_SRC/LOBSTER/trcsed.F90 (modified) (1 diff)
• TOP_SRC/PISCES/p4zche.F90 (modified) (3 diffs)
• TOP_SRC/PISCES/p4zflx.F90 (modified) (1 diff)
• TOP_SRC/PISCES/p4zint.F90 (modified) (1 diff)
• TOP_SRC/PISCES/p4zlim.F90 (modified) (1 diff)
• TOP_SRC/PISCES/p4zprod.F90 (modified) (1 diff)
• TOP_SRC/PISCES/p4zrem.F90 (modified) (1 diff)
• TOP_SRC/PISCES/p4zsed.F90 (modified) (3 diffs)
• TOP_SRC/SED/par_sed.F90 (modified) (1 diff)
• TOP_SRC/SED/sed.F90 (modified) (2 diffs)
• TOP_SRC/SED/seddta.F90 (modified) (1 diff)
• TOP_SRC/TRP/trcadv.F90 (modified) (3 diffs)
• TOP_SRC/TRP/trcldf.F90 (modified) (1 diff)
• TOP_SRC/TRP/trcnxt.F90 (modified) (3 diffs)
• TOP_SRC/TRP/trctrp.F90 (modified) (2 diffs)
• TOP_SRC/TRP/trczdf.F90 (modified) (1 diff)
• TOP_SRC/oce_trc.F90 (modified) (1 diff)
• TOP_SRC/trc.F90 (modified) (1 diff)
• TOP_SRC/trcini.F90 (modified) (6 diffs)

branches/DEV_r2006_merge_TRA_TRC/NEMO/NST_SRC/agrif_user.F90

@@ -307 +307 @@
 #if defined key_top
         ! Check passive tracer cell
-        IF( ndttrc .ne. 1 ) THEN
-           WRITE(*,*) 'ndttrc should be equal to 1'
+        IF( nn_dttrc .ne. 1 ) THEN
+           WRITE(*,*) 'nn_dttrc should be equal to 1'
         ENDIF
 #endif
@@ -444 +444 @@
          ENDIF
         ! Check passive tracer cell
-        IF( ndttrc .ne. 1 ) THEN
-           WRITE(*,*) 'ndttrc should be equal to 1'
+        IF( nn_dttrc .ne. 1 ) THEN
+           WRITE(*,*) 'nn_dttrc should be equal to 1'
         ENDIF
 

branches/DEV_r2006_merge_TRA_TRC/NEMO/OFF_SRC/dtadyn.F90

@@ -184 +184 @@
          
          IF( lk_ldfslp ) THEN
-            ! Computes slopes
-            ! Caution : here tn, sn and avt are used as workspace
-            tn (:,:,:) = tdta  (:,:,:,2)
-            sn (:,:,:) = sdta  (:,:,:,2)
-            avt(:,:,:) = avtdta(:,:,:,2)
-         
-            CALL eos( tn, sn, rhd, rhop )   ! Time-filtered in situ density 
-            CALL bn2( tn, sn, rn2 )         ! before Brunt-Vaisala frequency
+            ! Computes slopes. Caution : here tsn and avt are used as workspace
+            tsn (:,:,:,jp_tem) = tdta  (:,:,:,2)
+            tsn (:,:,:,jp_sal) = sdta  (:,:,:,2)
+            avt(:,:,:)         = avtdta(:,:,:,2)
+         
+            CALL eos( tsn, rhd, rhop )   ! Time-filtered in situ density 
+            CALL bn2( tsn, rn2 )         ! before Brunt-Vaisala frequency
             IF( ln_zps )   &
-               &   CALL zps_hde( kt, tn , sn , rhd,  &  ! Partial steps: before Horizontal DErivative
-               &                 gtu, gsu, gru,  &  ! of t, s, rd at the bottom ocean level
-               &                 gtv, gsv, grv )
+               &   CALL zps_hde( kt, jpts, tsn, gtsu, gtsv,  &  ! Partial steps: before Horizontal DErivative
+               &                           rhd, gru , grv   )    ! of t, s, rd at the bottom ocean level
                    CALL zdf_mxl( kt )              ! mixed layer depth
                    CALL ldf_slp( kt, rhd, rn2 )
@@ -213 +211 @@
          
          IF( lk_ldfslp ) THEN
-            ! Computes slopes
-            ! Caution : here tn, sn and avt are used as workspace
-            tn (:,:,:) = tdta  (:,:,:,2)
-            sn (:,:,:) = sdta  (:,:,:,2)
-            avt(:,:,:) = avtdta(:,:,:,2)
-
-            CALL eos( tn, sn, rhd, rhop )   ! Time-filtered in situ density 
-            CALL bn2( tn, sn, rn2 )         ! before Brunt-Vaisala frequency
+            ! Computes slopes. Caution : here tsn and avt are used as workspace
+            tsn (:,:,:,jp_tem) = tdta  (:,:,:,2)
+            tsn (:,:,:,jp_sal) = sdta  (:,:,:,2)
+            avt(:,:,:)         = avtdta(:,:,:,2)
+         
+            CALL eos( tsn, rhd, rhop )   ! Time-filtered in situ density 
+            CALL bn2( tsn, rn2 )         ! before Brunt-Vaisala frequency
             IF( ln_zps )   &
-               &   CALL zps_hde( kt, tn , sn , rhd,  &  ! Partial steps: before Horizontal DErivative
-               &                 gtu, gsu, gru,  &  ! of t, s, rd at the bottom ocean level
-               &                 gtv, gsv, grv )
+               &   CALL zps_hde( kt, jpts, tsn, gtsu, gtsv,  &  ! Partial steps: before Horizontal DErivative
+               &                           rhd, gru , grv   )    ! of t, s, rd at the bottom ocean level
                    CALL zdf_mxl( kt )              ! mixed layer depth
                    CALL ldf_slp( kt, rhd, rn2 )
@@ -261 +257 @@
 
          IF( lk_ldfslp ) THEN
-            ! Computes slopes
-            ! Caution : here tn, sn and avt are used as workspace
-            tn (:,:,:) = tdta  (:,:,:,2)
-            sn (:,:,:) = sdta  (:,:,:,2)
-            avt(:,:,:) = avtdta(:,:,:,2)
-
-            CALL eos( tn, sn, rhd, rhop )   ! Time-filtered in situ density 
-            CALL bn2( tn, sn, rn2 )         ! before Brunt-Vaisala frequency
+            ! Computes slopes. Caution : here tsn and avt are used as workspace
+            tsn (:,:,:,jp_tem) = tdta  (:,:,:,2)
+            tsn (:,:,:,jp_sal) = sdta  (:,:,:,2)
+            avt(:,:,:)         = avtdta(:,:,:,2)
+         
+            CALL eos( tsn, rhd, rhop )   ! Time-filtered in situ density 
+            CALL bn2( tsn, rn2 )         ! before Brunt-Vaisala frequency
             IF( ln_zps )   &
-               &   CALL zps_hde( kt, tn , sn , rhd,  &  ! Partial steps: before Horizontal DErivative
-               &                 gtu, gsu, gru,  &  ! of t, s, rd at the bottom ocean level
-               &                 gtv, gsv, grv )
+               &   CALL zps_hde( kt, jpts, tsn, gtsu, gtsv,  &  ! Partial steps: before Horizontal DErivative
+               &                           rhd, gru , grv   )    ! of t, s, rd at the bottom ocean level
                    CALL zdf_mxl( kt )              ! mixed layer depth
                    CALL ldf_slp( kt, rhd, rn2 )
@@ -311 +305 @@
       
       ! In any case, we need rhop
-      CALL eos( tn, sn, rhd, rhop ) 
+      CALL eos( tsn, rhd, rhop ) 
       
 #if ! defined key_degrad && defined key_traldf_c2d
@@ -320 +314 @@
       ! Compute bbl coefficients if needed
       IF( lk_trabbl ) THEN
-         tb(:,:,:) = tn(:,:,:)
-         sb(:,:,:) = sn(:,:,:)
+         tsb(:,:,:,:) = tsn(:,:,:,:)
          CALL bbl( kt, 'TRC')
       END IF
@@ -735 +728 @@
       !!----------------------------------------------------------------------
       
-      tn (:,:,:) = tdta  (:,:,:,2)
-      sn (:,:,:) = sdta  (:,:,:,2)
-      avt(:,:,:) = avtdta(:,:,:,2)
+      tsn(:,:,:,jp_tem) = tdta  (:,:,:,2)
+      tsn(:,:,:,jp_sal) = sdta  (:,:,:,2)
+      avt(:,:,:)        = avtdta(:,:,:,2)
       
       un (:,:,:) = udta  (:,:,:,2) 
@@ -796 +789 @@
       zweighm1 = 1. - pweigh
       
-      tn (:,:,:) = zweighm1 * tdta  (:,:,:,1) + pweigh * tdta  (:,:,:,2)
-      sn (:,:,:) = zweighm1 * sdta  (:,:,:,1) + pweigh * sdta  (:,:,:,2)
-      avt(:,:,:) = zweighm1 * avtdta(:,:,:,1) + pweigh * avtdta(:,:,:,2)
+      tsn(:,:,:,jp_tem) = zweighm1 * tdta  (:,:,:,1) + pweigh * tdta  (:,:,:,2)
+      tsn(:,:,:,jp_sal) = zweighm1 * sdta  (:,:,:,1) + pweigh * sdta  (:,:,:,2)
+      avt(:,:,:)        = zweighm1 * avtdta(:,:,:,1) + pweigh * avtdta(:,:,:,2)
       
       un (:,:,:) = zweighm1 * udta  (:,:,:,1) + pweigh * udta  (:,:,:,2) 

branches/DEV_r2006_merge_TRA_TRC/NEMO/OFF_SRC/istate.F90

@@ -59 +59 @@
       ;   hdivn(:,:,:) = 0.e0   ;
 
-      ;   tn   (:,:,:) = 0.e0   ;   ta   (:,:,:) = 0.e0
-      ;   sn   (:,:,:) = 0.e0   ;   sa   (:,:,:) = 0.e0
+      ;   tsn  (:,:,:,:) = 0.e0 
 
       rhd  (:,:,:) = 0.e0

branches/DEV_r2006_merge_TRA_TRC/NEMO/OFF_SRC/opa.F90

@@ -223 +223 @@
 
 
+      !                                     ! General initialization
                             CALL     phy_cst    ! Physical constants
                             CALL     eos_init   ! Equation of state
@@ -234 +235 @@
       IF( lk_zdfddm .AND. .NOT. lk_zdfkpp )   &
          &                  CALL zdf_ddm_init   ! double diffusive mixing
-      !                                     ! Lateral physics
 #if ! defined key_degrad
                             CALL ldf_tra_init   ! Lateral ocean tracer physics
 #endif
       IF( lk_ldfslp )       CALL ldf_slp_init   ! slope of lateral mixing
+
       !                                     ! Active tracers
                             CALL tra_qsr_init   ! penetrative solar radiation qsr
       IF( lk_trabbl     )   CALL tra_bbl_init   ! advective (and/or diffusive) bottom boundary layer scheme
 
-                            CALL     trc_ini    ! Passive tracers
+      !                                     ! Passive tracers
+                            CALL     trc_init   ! Passive tracers initialization
+      !                                     ! Dynamics
                             CALL dta_dyn_init   ! Initialization for the dynamics
                             CALL     iom_init       ! iom_put initialization

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/DIA/diaar5.F90

@@ -67 +67 @@
       REAL(wp), DIMENSION(jpi,jpj    ) ::   zarea_ssh, zbotpres
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zrhd, zrhop
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) ::   ztsn
       !!--------------------------------------------------------------------
 
@@ -82 +83 @@
 
       !                                         ! thermosteric ssh
-      CALL eos( tn, sn0, zrhd )                       ! now in situ density using initial salinity
+      ztsn(:,:,:,jp_tem) = tn (:,:,:)
+      ztsn(:,:,:,jp_sal) = sn0(:,:,:)
+      CALL eos( ztsn, zrhd )                       ! now in situ density using initial salinity
       !
       zbotpres(:,:) = 0.e0                            ! no atmospheric surface pressure, levitating sea-ice
@@ -96 +99 @@
       
       !                                         ! steric sea surface height
-      CALL eos( tn, sn, zrhd, zrhop )                 ! now in situ and potential density
+      CALL eos( tsn, zrhd, zrhop )                 ! now in situ and potential density
       zrhop(:,:,jpk) = 0.e0
       CALL iom_put( 'rhop', zrhop )

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/IOM/restart.F90

@@ -26 +26 @@
    USE zdfmxl          ! mixed layer depth
    USE trdmld_oce      ! ocean active mixed layer tracers trends variables
-
+#if defined key_zdfkpp
+   USE traswap
+#endif
    IMPLICIT NONE
    PRIVATE
@@ -205 +207 @@
          CALL iom_get( numror, jpdom_autoglo, 'rhd' , rhd  )       ! now    in situ density anomaly
       ELSE
-         CALL eos( tn, sn, rhd )   ! compute rhd
+         CALL tra_swap
+         CALL eos( tsn, rhd )   ! compute rhd
       ENDIF
 #endif

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/eosbn2.F90

@@ -68 +68 @@
 CONTAINS
 
-   SUBROUTINE eos_insitu( ptem, psal, prd )
+   SUBROUTINE eos_insitu( pts, prd )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE eos_insitu  ***
@@ -98 +98 @@
       !!              prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
       !!      Note that no boundary condition problem occurs in this routine
-      !!      as (ptem,psal) are defined over the whole domain.
+      !!      as pts are defined over the whole domain.
       !!
       !! ** Action  :   compute prd , the in situ density (no units)
@@ -104 +104 @@
       !! References :   Jackett and McDougall, J. Atmos. Ocean. Tech., 1994
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   ptem   ! potential temperature  [Celcius]
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   psal   ! salinity               [psu]
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   prd    ! in situ density 
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius]
+      !                                                               ! 2 : salinity               [psu]
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   prd   ! in situ density 
       !!
       INTEGER  ::   ji, jj, jk           ! dummy loop indices
@@ -123 +123 @@
          zrau0r = 1.e0 / rau0
 !CDIR NOVERRCHK
-         zws(:,:,:) = SQRT( ABS( psal(:,:,:) ) )
+         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
          !  
          DO jk = 1, jpkm1
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zt = ptem  (ji,jj,jk)
-                  zs = psal  (ji,jj,jk)
+                  zt = pts   (ji,jj,jk,jp_tem)
+                  zs = pts   (ji,jj,jk,jp_sal)
                   zh = fsdept(ji,jj,jk)        ! depth
                   zsr= zws   (ji,jj,jk)        ! square root salinity
@@ -169 +169 @@
       CASE( 1 )                !==  Linear formulation function of temperature only  ==!
          DO jk = 1, jpkm1
-            prd(:,:,jk) = ( 0.0285 - rn_alpha * ptem(:,:,jk) ) * tmask(:,:,jk)
+            prd(:,:,jk) = ( 0.0285 - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk)
          END DO
          !
       CASE( 2 )                !==  Linear formulation function of temperature and salinity  ==!
          DO jk = 1, jpkm1
-            prd(:,:,jk) = ( rn_beta  * psal(:,:,jk) - rn_alpha * ptem(:,:,jk) ) * tmask(:,:,jk)
+            prd(:,:,jk) = ( rn_beta  * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk)
          END DO
          !
@@ -184 +184 @@
 
 
-   SUBROUTINE eos_insitu_pot( ptem, psal, prd, prhop )
+   SUBROUTINE eos_insitu_pot( pts, prd, prhop )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE eos_insitu_pot  ***
@@ -230 +230 @@
       !!                Brown and Campana, Mon. Weather Rev., 1978
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   ptem   ! potential temperature  [Celcius]
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   psal   ! salinity               [psu]
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   prd    ! in situ density 
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   prhop  ! potential density (surface referenced)
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius]
+      !                                                               ! 2 : salinity               [psu]
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd    ! in situ density 
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prhop  ! potential density (surface referenced)
 
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -246 +246 @@
          zrau0r = 1.e0 / rau0
 !CDIR NOVERRCHK
-         zws(:,:,:) = SQRT( ABS( psal(:,:,:) ) )
+         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
          !  
          DO jk = 1, jpkm1
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zt = ptem  (ji,jj,jk)
-                  zs = psal  (ji,jj,jk)
+                  zt = pts   (ji,jj,jk,jp_tem)
+                  zs = pts   (ji,jj,jk,jp_sal)
                   zh = fsdept(ji,jj,jk)        ! depth
                   zsr= zws   (ji,jj,jk)        ! square root salinity
@@ -295 +295 @@
       CASE( 1 )                !==  Linear formulation = F( temperature )  ==!
          DO jk = 1, jpkm1
-            prd  (:,:,jk) = ( 0.0285 - rn_alpha * ptem(:,:,jk) )        * tmask(:,:,jk)
-            prhop(:,:,jk) = ( 1.e0   +          prd (:,:,jk) ) * rau0 * tmask(:,:,jk)
+            prd  (:,:,jk) = ( 0.0285 - rn_alpha * pts(:,:,jk,jp_sal) )        * tmask(:,:,jk)
+            prhop(:,:,jk) = ( 1.e0   +            prd (:,:,jk)       ) * rau0 * tmask(:,:,jk)
          END DO
          !
       CASE( 2 )                !==  Linear formulation = F( temperature , salinity )  ==!
          DO jk = 1, jpkm1
-            prd  (:,:,jk) = ( rn_beta  * psal(:,:,jk) - rn_alpha * ptem(:,:,jk) )        * tmask(:,:,jk)
-            prhop(:,:,jk) = ( 1.e0   + prd (:,:,jk) )                         * rau0 * tmask(:,:,jk)
+            prd  (:,:,jk) = ( rn_beta  * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) )        * tmask(:,:,jk)
+            prhop(:,:,jk) = ( 1.e0     + prd (:,:,jk) )                                       * rau0 * tmask(:,:,jk)
          END DO
          !
@@ -312 +312 @@
 
 
-   SUBROUTINE eos_insitu_2d( ptem, psal, pdep, prd )
+   SUBROUTINE eos_insitu_2d( pts, pdep, prd )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE eos_insitu_2d  ***
@@ -342 +342 @@
       !!              prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
       !!      Note that no boundary condition problem occurs in this routine
-      !!      as (ptem,psal) are defined over the whole domain.
+      !!      as pts are defined over the whole domain.
       !!
       !! ** Action  : - prd , the in situ density (no units)
@@ -348 +348 @@
       !! References :   Jackett and McDougall, J. Atmos. Ocean. Tech., 1994
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptem   ! potential temperature  [Celcius]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity               [psu]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pdep   ! depth                  [m]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   prd    ! in situ density 
+      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius]
+      !                                                           ! 2 : salinity               [psu]
+      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::   pdep  ! depth                  [m]
+      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density 
       !!
       INTEGER  ::   ji, jj                    ! dummy loop indices
@@ -369 +369 @@
 !CDIR NOVERRCHK
             DO ji = 1, fs_jpim1   ! vector opt.
-               zws(ji,jj) = SQRT( ABS( psal(ji,jj) ) )
+               zws(ji,jj) = SQRT( ABS( pts(ji,jj,jp_sal) ) )
             END DO
          END DO
@@ -375 +375 @@
             DO ji = 1, fs_jpim1   ! vector opt.
                zmask = tmask(ji,jj,1)          ! land/sea bottom mask = surf. mask
-               zt    = ptem (ji,jj)            ! interpolated T
-               zs    = psal (ji,jj)            ! interpolated S
+               zt    = pts  (ji,jj,jp_tem)            ! interpolated T
+               zs    = pts  (ji,jj,jp_sal)            ! interpolated S
                zsr   = zws  (ji,jj)            ! square root of interpolated S
                zh    = pdep (ji,jj)            ! depth at the partial step level
@@ -416 +416 @@
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               prd(ji,jj) = ( 0.0285 - rn_alpha * ptem(ji,jj) ) * tmask(ji,jj,1)
+               prd(ji,jj) = ( 0.0285 - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1)
             END DO
          END DO
@@ -423 +423 @@
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               prd(ji,jj) = ( rn_beta * psal(ji,jj) - rn_alpha * ptem(ji,jj) ) * tmask(ji,jj,1) 
+               prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 
             END DO
          END DO
@@ -434 +434 @@
 
 
-   SUBROUTINE eos_bn2( ptem, psal, pn2 )
+   SUBROUTINE eos_bn2( pts, pn2 )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE eos_bn2  ***
@@ -464 +464 @@
       !! References :   McDougall, J. Phys. Oceanogr., 17, 1950-1964, 1987.
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   ptem   ! potential temperature   [Celcius]
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   psal   ! salinity                [psu]
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   pn2    ! Brunt-Vaisala frequency [s-1]
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius]
+      !                                                               ! 2 : salinity               [psu]
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pn2    ! Brunt-Vaisala frequency [s-1]
       !!
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -485 +485 @@
                DO ji = 1, jpi
                   zgde3w = grav / fse3w(ji,jj,jk)
-                  zt = 0.5 * ( ptem(ji,jj,jk) + ptem(ji,jj,jk-1) )          ! potential temperature at w-point
-                  zs = 0.5 * ( psal(ji,jj,jk) + psal(ji,jj,jk-1) ) - 35.0   ! salinity anomaly (s-35) at w-point
+                  zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) )          ! potential temperature at w-point
+                  zs = 0.5 * ( pts(ji,jj,jk,jp_sal) + pts(ji,jj,jk-1,jp_sal) ) - 35.0   ! salinity anomaly (s-35) at w-point
                   zh = fsdepw(ji,jj,jk)                                     ! depth in meters  at w-point
                   !
@@ -515 +515 @@
                      !
                   pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk)           &   ! N^2 
-                     &          * ( zalbet * ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) )   &
-                     &                     - ( psal(ji,jj,jk-1) - psal(ji,jj,jk) ) )
+                     &          * ( zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   &
+                     &                     - ( pts(ji,jj,jk-1njp_sal) - pts(ji,jj,jk,jp_sal) ) )
 #if defined key_zdfddm
                   !                                                         !!bug **** caution a traiter zds=dk[S]= 0 !!!!
-                  zds = ( psal(ji,jj,jk-1) - psal(ji,jj,jk) )                    ! Rrau = (alpha / beta) (dk[t] / dk[s])
+                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )                    ! Rrau = (alpha / beta) (dk[t] / dk[s])
                   IF ( ABS( zds) <= 1.e-20 ) zds = 1.e-20
-                  rrau(ji,jj,jk) = zalbet * ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) / zds
+                  rrau(ji,jj,jk) = zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
 #endif
                END DO
@@ -529 +529 @@
       CASE( 1 )                !==  Linear formulation = F( temperature )  ==!
          DO jk = 2, jpkm1
-            pn2(:,:,jk) = grav * rn_alpha * ( ptem(:,:,jk-1) - ptem(:,:,jk) ) / fse3w(:,:,jk) * tmask(:,:,jk)
+            pn2(:,:,jk) = grav * rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) ) / fse3w(:,:,jk) * tmask(:,:,jk)
          END DO
          !
       CASE( 2 )                !==  Linear formulation = F( temperature , salinity )  ==!
          DO jk = 2, jpkm1
-            pn2(:,:,jk) = grav * (  rn_alpha * ( ptem(:,:,jk-1) - ptem(:,:,jk) )      &
-               &                  - rn_beta  * ( psal(:,:,jk-1) - psal(:,:,jk) )  )   &
+            pn2(:,:,jk) = grav * (  rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) )      &
+               &                  - rn_beta  * ( pts(:,:,jk-1,jp_sal) - pts(:,:,jk,jp_sal) )  )   &
                &               / fse3w(:,:,jk) * tmask(:,:,jk)
          END DO 
@@ -542 +542 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zds = ( psal(ji,jj,jk-1) - psal(ji,jj,jk) )
+                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  
                   IF ( ABS( zds ) <= 1.e-20 ) zds = 1.e-20
-                  rrau(ji,jj,jk) = ralpbet * ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) / zds
+                  rrau(ji,jj,jk) = ralpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
                END DO
             END DO

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv.F90

@@ -42 +42 @@
    INTEGER ::   nadv   ! choice of the type of advection scheme
 
+   REAL(wp), DIMENSION(jpk) ::   r2dt  ! vertical profile time-step, = 2 rdttra
+      !                                ! except at nit000 (=rdttra) if neuler=0
+
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -66 +69 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk)   ::  zun, zvn, zwn   ! effective transport
       !!----------------------------------------------------------------------
+      !                                          ! set time step
+      IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000
+         r2dt(:) =  rdttra(:)                          ! = rdtra (restarting with Euler time stepping)
+      ELSEIF( kt <= nit000 + 1) THEN                ! at nit000 or nit000+1
+         r2dt(:) = 2. * rdttra(:)                      ! = 2 rdttra (leapfrog)
+      ENDIF
       !                                               !==  effective transport  ==!
       DO jk = 1, jpkm1
@@ -82 +91 @@
 
       SELECT CASE ( nadv )                            !==  compute advection trend and add it to general trend  ==!
-      CASE ( 1 )   ;    CALL tra_adv_cen2  ( kt, 'TRA', zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  2nd order centered
-      CASE ( 2 )   ;    CALL tra_adv_tvd   ( kt, 'TRA', zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  TVD 
-      CASE ( 3 )   ;    CALL tra_adv_muscl ( kt, 'TRA', zun, zvn, zwn, tsb,      tsa, jpts )   !  MUSCL 
-      CASE ( 4 )   ;    CALL tra_adv_muscl2( kt, 'TRA', zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  MUSCL2 
-      CASE ( 5 )   ;    CALL tra_adv_ubs   ( kt, 'TRA', zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  UBS 
-      CASE ( 6 )   ;    CALL tra_adv_qck   ( kt, 'TRA', zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  QUICKEST 
+      CASE ( 1 )   ;    CALL tra_adv_cen2  ( kt, 'TRA',       zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  2nd order centered
+      CASE ( 2 )   ;    CALL tra_adv_tvd   ( kt, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  TVD 
+      CASE ( 3 )   ;    CALL tra_adv_muscl ( kt, 'TRA', r2dt, zun, zvn, zwn, tsb,      tsa, jpts )   !  MUSCL 
+      CASE ( 4 )   ;    CALL tra_adv_muscl2( kt, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  MUSCL2 
+      CASE ( 5 )   ;    CALL tra_adv_ubs   ( kt, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  UBS 
+      CASE ( 6 )   ;    CALL tra_adv_qck   ( kt, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  QUICKEST 
       !
       CASE (-1 )                                      !==  esopa: test all possibility with control print  ==!
-         CALL tra_adv_cen2  ( kt , 'TRA', zun, zvn, zwn, tsb, tsn, tsa, jpts )          
+         CALL tra_adv_cen2  ( kt, 'TRA',       zun, zvn, zwn, tsb, tsn, tsa, jpts )          
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv0 - Ta: ', mask1=tmask,               &
             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_adv_tvd   ( kt, 'TRA', zun, zvn, zwn, tsb, tsn, tsa, jpts )          
+         CALL tra_adv_tvd   ( kt, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts )          
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv1 - Ta: ', mask1=tmask,               &
             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_adv_muscl ( kt, 'TRA', zun, zvn, zwn, tsb,      tsa, jpts )          
+         CALL tra_adv_muscl ( kt, 'TRA', r2dt, zun, zvn, zwn, tsb,      tsa, jpts )          
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv3 - Ta: ', mask1=tmask,               &
             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_adv_muscl2( kt, 'TRA', zun, zvn, zwn, tsb, tsn, tsa, jpts )          
+         CALL tra_adv_muscl2( kt, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts )          
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv4 - Ta: ', mask1=tmask,               &
             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_adv_ubs   ( kt, 'TRA', zun, zvn, zwn, tsb, tsn, tsa, jpts )          
+         CALL tra_adv_ubs   ( kt, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts )          
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv5 - Ta: ', mask1=tmask,               &
             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_adv_qck   ( kt, 'TRA', zun, zvn, zwn, tsb, tsn, tsa, jpts )          
+         CALL tra_adv_qck   ( kt, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts )          
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv6 - Ta: ', mask1=tmask,               &
             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
@@ -179 +188 @@
          IF( nadv ==  5 )   WRITE(numout,*) '         UBS       scheme is used'
          IF( nadv ==  6 )   WRITE(numout,*) '         QUICKEST  scheme is used'
-         IF( nadv ==  7 )   WRITE(numout,*) '         SMOLAR    scheme is used'
          IF( nadv == -1 )   WRITE(numout,*) '         esopa test: use all advection scheme'
       ENDIF

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_cen2.F90

@@ -20 +20 @@
    !!                  area (set for orca 2 and 4 only)
    !!----------------------------------------------------------------------
-   USE oce, ONLY: tn, sn  ! now ocean temperature and salinity
+   USE oce, ONLY: tsn  ! now ocean temperature and salinity
    USE dom_oce         ! ocean space and time domain
    USE eosbn2          ! equation of state
@@ -32 +32 @@
    USE zdf_oce         ! ocean vertical physics
    USE restart         ! ocean restart
+   USE trc_oce         ! share passive tracers/Ocean variables
 
    IMPLICIT NONE
@@ -134 +135 @@
       !!----------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
+
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme'
+         IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme on ', cdtype
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~   Vector optimization case'
          IF(lwp) WRITE(numout,*)
@@ -159 +161 @@
 !!gm  not strickly exact : the freezing point should be computed at each ocean levels...
 !!gm  not a big deal since cen2 is no more used in global ice-ocean simulations
-      ztfreez(:,:) = tfreez( sn(:,:,1) )
+      ztfreez(:,:) = tfreez( tsn(:,:,1,jp_sal) )
       DO jk = 1, jpk
          DO jj = 1, jpj
             DO ji = 1, jpi
                !                                        ! below ice covered area (if tn < "freezing"+0.1 )
-               IF( tn(ji,jj,jk) <= ztfreez(ji,jj) + 0.1 ) THEN   ;   zice = 1.e0
-               ELSE                                              ;   zice = 0.e0
+               IF( tsn(ji,jj,jk,jp_tem) <= ztfreez(ji,jj) + 0.1 ) THEN   ;   zice = 1.e0
+               ELSE                                                      ;   zice = 0.e0
                ENDIF
                zind(ji,jj,jk) = MAX (   &

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_eiv.F90

@@ -22 +22 @@
    USE in_out_manager  ! I/O manager
    USE iom
+   USE trc_oce         ! share passive tracers/Ocean variables
 # if defined key_diaeiv
    USE phycst          ! physical constants
@@ -80 +81 @@
       !!----------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_adv_eiv : eddy induced advection :'
+         IF(lwp) WRITE(numout,*) 'tra_adv_eiv : eddy induced advection on ', cdtype,' :'
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   add to velocity fields the eiv component'
 # if defined key_diaeiv 

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_muscl.F90

@@ -24 +24 @@
    USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    USE diaptr          ! poleward transport diagnostics
+   USE trc_oce         ! share passive tracers/Ocean variables
 
 
@@ -44 +45 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_muscl( kt, cdtype, pun, pvn, pwn, &
-      &                                       ptb, pta, kjpt   )
+   SUBROUTINE tra_adv_muscl( kt, cdtype, p2dt, pun, pvn, pwn, &
+      &                                        ptb, pta, kjpt   )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE tra_adv_muscl  ***
@@ -67 +68 @@
       CHARACTER(len=3), INTENT(in   )                               ::   cdtype          ! =TRA or TRC (tracer indicator)
       INTEGER         , INTENT(in   )                               ::   kjpt            ! number of tracers
+      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::   p2dt            ! vertical profile of tracer time-step
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pun, pvn, pwn   ! 3 ocean velocity components
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb           ! before and now tracer fields
@@ -75 +77 @@
       REAL(wp) ::   zv, z0v, zzwy
       REAL(wp) ::   zw, z0w
-      REAL(wp) ::   ztra, zbtr, z2, zdt, zalpha
+      REAL(wp) ::   ztra, zbtr, zdt, zalpha
       REAL(wp), DIMENSION (jpi,jpj,jpk) :: zslpx, zslpy   ! 3D workspace
       !!----------------------------------------------------------------------
 
-      IF( kt == nit000 .AND. lwp ) THEN
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          WRITE(numout,*)
-         WRITE(numout,*) 'tra_adv : MUSCL advection scheme'
+         WRITE(numout,*) 'tra_adv : MUSCL advection scheme on ', cdtype
          WRITE(numout,*) '~~~~~~~'
          !
@@ -88 +90 @@
       ENDIF
 
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.
-      ELSE                                        ;    z2 = 2.
-      ENDIF
-      !
       !                                                     ! ===========
       DO jn = 1, kjpt                                       ! tracer loop
@@ -139 +137 @@
          !                                             !-- MUSCL horizontal advective fluxes
          DO jk = 1, jpkm1                                     ! interior values
-            zdt  = z2 * rdttra(jk)
+            zdt  = p2dt(jk)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -230 +228 @@
          !
          DO jk = 1, jpkm1                                     ! interior values
-            zdt  = z2 * rdttra(jk)
+            zdt  = p2dt(jk)
             DO jj = 2, jpjm1      
                DO ji = fs_2, fs_jpim1   ! vector opt.

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_muscl2.F90

@@ -22 +22 @@
    USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    USE diaptr          ! poleward transport diagnostics
+   USE trc_oce         ! share passive tracers/Ocean variables
 
 
@@ -43 +44 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_muscl2( kt, cdtype, pun, pvn, pwn, &
-      &                                   ptb, ptn, pta, kjpt   )
+   SUBROUTINE tra_adv_muscl2( kt, cdtype, p2dt, pun, pvn, pwn, &
+      &                                         ptb, ptn, pta, kjpt   )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE tra_adv_muscl2  ***
@@ -67 +68 @@
       CHARACTER(len=3), INTENT(in   )                               ::   cdtype          ! =TRA or TRC (tracer indicator)
       INTEGER         , INTENT(in   )                               ::   kjpt            ! number of tracers
+      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::   p2dt            ! vertical profile of tracer time-step
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pun, pvn, pwn   ! 3 ocean velocity components
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb, ptn        ! before and now tracer fields
@@ -75 +77 @@
       REAL(wp) ::   zv, z0v, zzwy
       REAL(wp) ::   zw, z0w
-      REAL(wp) ::   ztra, zbtr, z2, zdt, zalpha
+      REAL(wp) ::   ztra, zbtr, zdt, zalpha
       REAL(wp), DIMENSION (jpi,jpj,jpk) ::  zslpx, zslpy   ! 3D workspace
       !!----------------------------------------------------------------------
 
-      IF( kt == nit000 .AND. lwp ) THEN
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          WRITE(numout,*)
-         WRITE(numout,*) 'tra_adv_muscl2 : MUSCL2 advection scheme'
+         WRITE(numout,*) 'tra_adv_muscl2 : MUSCL2 advection scheme on ', cdtype
          WRITE(numout,*) '~~~~~~~~~~~~~~~'
          !
@@ -88 +90 @@
       ENDIF
 
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   z2 = 1.
-      ELSE                                        ;   z2 = 2.
-      ENDIF
-      !
       !
       DO jn = 1, kjpt                                            ! tracer loop
@@ -139 +137 @@
         !                                             !-- MUSCL horizontal advective fluxes
          DO jk = 1, jpkm1                                     ! interior values
-            zdt  = z2 * rdttra(jk)
+            zdt  = p2dt(jk)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -255 +253 @@
          !
          DO jk = 1, jpkm1                                     ! interior values
-            zdt  = z2 * rdttra(jk)
+            zdt  = p2dt(jk)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_qck.F90

@@ -25 +25 @@
    USE in_out_manager  ! I/O manager
    USE diaptr          ! poleward transport diagnostics
+   USE trc_oce         ! share passive tracers/Ocean variables
 
    IMPLICIT NONE
@@ -45 +46 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_qck ( kt, cdtype, pun, pvn, pwn, &
-      &                                 ptb, ptn, pta, kjpt   )
+   SUBROUTINE tra_adv_qck ( kt, cdtype, p2dt, pun, pvn, pwn, &
+      &                                       ptb, ptn, pta, kjpt   )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv_qck  ***
@@ -85 +86 @@
       CHARACTER(len=3), INTENT(in   )                               ::   cdtype          ! =TRA or TRC (tracer indicator)
       INTEGER         , INTENT(in   )                               ::   kjpt            ! number of tracers
+      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::   p2dt            ! vertical profile of tracer time-step
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pun, pvn, pwn   ! 3 ocean velocity components
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb, ptn        ! before and now tracer fields
       REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta           ! tracer trend 
-      !!
-      REAL(wp) ::   z2                            !    temporary scalar
-      !!----------------------------------------------------------------------
-
-      IF( kt == nit000 ) THEN
+      !!----------------------------------------------------------------------
+
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme'
+         IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
          IF(lwp) WRITE(numout,*)
@@ -102 +102 @@
       ENDIF
 
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.
-      ELSE                                        ;    z2 = 2.
-      ENDIF
-
       ! I. The horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
       !---------------------------------------------------------------------------
 
-      CALL tra_adv_qck_i( kt, cdtype, pun, z2, ptb, ptn, pta, kjpt ) 
-      CALL tra_adv_qck_j( kt, cdtype, pvn, z2, ptb, ptn, pta, kjpt ) 
+      CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt ) 
+      CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt ) 
 
       ! II. The vertical fluxes are computed with the 2nd order centered scheme
@@ -119 +115 @@
    END SUBROUTINE tra_adv_qck
 
-   SUBROUTINE tra_adv_qck_i( kt, cdtype, pun, pz2,   &
-      &                                  ptb, ptn, pta, kjpt   )
+   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun,    &
+      &                                        ptb, ptn, pta, kjpt   )
       !!----------------------------------------------------------------------
       !!
@@ -129 +125 @@
       CHARACTER(len=3), INTENT(in   )                               ::   cdtype          ! =TRA or TRC (tracer indicator)
       INTEGER         , INTENT(in   )                               ::   kjpt            ! number of tracers
+      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::   p2dt            ! vertical profile of tracer time-step
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pun             ! zonal velocity component
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb, ptn    ! before tracer fields
       REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta           ! tracer trend 
-      REAL(wp)        , INTENT(in   )                               ::   pz2
       !!
       INTEGER  :: ji, jj, jk, jn           ! dummy loop indices
@@ -176 +172 @@
          !
          DO jk = 1, jpkm1  
-            zdt =  pz2 * rdttra(jk)
+            zdt =  p2dt(jk)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.   
@@ -239 +235 @@
    END SUBROUTINE tra_adv_qck_i
 
-   SUBROUTINE tra_adv_qck_j( kt, cdtype, pvn, pz2,   &
-      &                                  ptb, ptn, pta, kjpt   )
+   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn,   &
+      &                                        ptb, ptn, pta, kjpt   )
       !!----------------------------------------------------------------------
       !!
@@ -250 +246 @@
       CHARACTER(len=3), INTENT(in   )                               ::   cdtype          ! =TRA or TRC (tracer indicator)
       INTEGER         , INTENT(in   )                               ::   kjpt            ! number of tracers
+      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::   p2dt            ! vertical profile of tracer time-step
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pvn             ! meridional velocity component
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb, ptn    ! before tracer fields
       REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta           ! tracer trend 
-      REAL(wp)        , INTENT(in   )                               ::   pz2
       !!
       INTEGER  :: ji, jj, jk, jn           ! dummy loop indices
@@ -296 +292 @@
          !
          DO jk = 1, jpkm1  
-            zdt =  pz2 * rdttra(jk)
+            zdt =  p2dt(jk)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.   

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_tvd.F90

@@ -32 +32 @@
    USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    USE diaptr          ! poleward transport diagnostics
+   USE trc_oce         ! share passive tracers/Ocean variables
 
 
@@ -52 +53 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_tvd ( kt, cdtype, pun, pvn, pwn, &
-      &                                 ptb, ptn, pta, kjpt   )
+   SUBROUTINE tra_adv_tvd ( kt, cdtype, p2dt, pun, pvn, pwn, &
+      &                                       ptb, ptn, pta, kjpt   )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv_tvd  ***
@@ -73 +74 @@
       CHARACTER(len=3), INTENT(in   )                               ::   cdtype          ! =TRA or TRC (tracer indicator)
       INTEGER         , INTENT(in   )                               ::   kjpt            ! number of tracers
+      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::   p2dt            ! vertical profile of tracer time-step
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pun, pvn, pwn   ! 3 ocean velocity components
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb, ptn        ! before and now tracer fields
@@ -79 +81 @@
       INTEGER  ::   ji, jj, jk, jn          ! dummy loop indices  
       REAL(wp) ::   &
-         z2, z2dtt, zbtr, ztra,              &  ! temporary scalar
+         z2dtt, zbtr, ztra,                  &  ! temporary scalar
          zfp_ui, zfp_vj, zfp_wk,             &  !    "         "
          zfm_ui, zfm_vj, zfm_wk                 !    "         "
@@ -86 +88 @@
       !!----------------------------------------------------------------------
 
-      zwi(:,:,:) = 0.e0
-
-      IF( kt == nit000 .AND. lwp ) THEN
+
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          WRITE(numout,*)
-         WRITE(numout,*) 'tra_adv_tvd : TVD advection scheme'
+         WRITE(numout,*) 'tra_adv_tvd : TVD advection scheme on ', cdtype
          WRITE(numout,*) '~~~~~~~~~~~'
          !
@@ -103 +104 @@
       END IF
       !
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.
-      ELSE                                        ;    z2 = 2.
-      ENDIF
+      zwi(:,:,:) = 0.e0
       !
       !                                                          ! ===========
@@ -150 +149 @@
          ! total advective trend
          DO jk = 1, jpkm1
-            z2dtt = z2 * rdttra(jk)
+            z2dtt = p2dt(jk)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -209 +208 @@
          ! 4. monotonicity algorithm
          ! -------------------------
-         CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, z2 )
+         CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
 
 
@@ -253 +252 @@
 
 
-   SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, prdt )
+   SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE nonosc  ***
@@ -266 +265 @@
       !!       in-space based differencing for fluid
       !!----------------------------------------------------------------------
-      REAL(wp), INTENT( in ) ::   prdt   ! ???
+      REAL(wp), DIMENSION(jpk)         , INTENT( in    ) ::   &  
+         p2dt                               ! vertical profile of tracer time-step
       REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT( in    ) ::   &
          pbef,                            & ! before field
@@ -299 +299 @@
       DO jk = 1, jpkm1
          ikm1 = MAX(jk-1,1)
-         z2dtt = prdt * rdttra(jk)
+         z2dtt = p2dt(jk)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_ubs.F90

@@ -21 +21 @@
    USE diaptr          ! poleward transport diagnostics
    USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
+   USE trc_oce         ! share passive tracers/Ocean variables
 
    IMPLICIT NONE
@@ -40 +41 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_ubs ( kt, cdtype, pun, pvn, pwn, &
-      &                                 ptb, ptn, pta, kjpt   )
+   SUBROUTINE tra_adv_ubs ( kt, cdtype, p2dt, pun, pvn, pwn, &
+      &                                       ptb, ptn, pta, kjpt   )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv_ubs  ***
@@ -81 +82 @@
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pun, pvn, pwn   ! 3 ocean velocity components
       INTEGER         , INTENT(in   )                               ::   kjpt            ! number of tracers
+      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::   p2dt            ! vertical profile of tracer time-step
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb, ptn        ! before and now tracer fields
       REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta           ! tracer trend 
@@ -88 +90 @@
       REAL(wp) ::   zfp_ui, zfm_ui, zcenut  !    "         "
       REAL(wp) ::   zfp_vj, zfm_vj, zcenvt  !    "         "    !    "         "
-      REAL(wp) ::   z2dtt, z2                   
+      REAL(wp) ::   z2dtt                   
       REAL(wp) ::   ztak, zfp_wk, zfm_wk    !    "         "
       REAL(wp) ::   zeeu, zeev, z_hdivn     
@@ -96 +98 @@
 
 
-      IF( kt == nit000 ) THEN
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_adv_ubs :  horizontal UBS advection scheme'
+         IF(lwp) WRITE(numout,*) 'tra_adv_ubs :  horizontal UBS advection scheme on ', cdtype
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
          !
          l_trd = .FALSE.
          IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
-      ENDIF
-      !
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.
-      ELSE                                        ;    z2 = 2.
       ENDIF
       !
@@ -232 +230 @@
          ! update and guess with monotonic sheme
          DO jk = 1, jpkm1
-            z2dtt = z2 * rdttra(jk)
+            z2dtt = p2dt(jk)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -255 +253 @@
          END DO
          !
-         CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, z2 )      !  monotonicity algorithm
+         CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt )      !  monotonicity algorithm
 
          !  final trend with corrected fluxes
@@ -288 +286 @@
    END SUBROUTINE tra_adv_ubs
 
-   SUBROUTINE nonosc_z( pbef, pcc, paft, prdt )
+   SUBROUTINE nonosc_z( pbef, pcc, paft, p2dt )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE nonosc_z  ***
@@ -301 +299 @@
       !!       in-space based differencing for fluid
       !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in   )                          ::   prdt   ! ???
+      REAL(wp), INTENT(in   ), DIMENSION(jpk)          ::   p2dt            ! vertical profile of tracer time-step
       REAL(wp),                DIMENSION (jpi,jpj,jpk) ::   pbef   ! before field
       REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   paft   ! after field
@@ -356 +354 @@
 
       DO jk = 1, jpkm1
-         z2dtt = prdt * rdttra(jk)
+         z2dtt = p2dt(jk)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -374 +372 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-
                za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
                zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trabbl.F90

@@ -32 +32 @@
    USE lbclnk         ! ocean lateral boundary conditions
    USE prtctl         ! Print control
+   USE trc_oce         ! share passive tracers/Ocean variables
 
    IMPLICIT NONE
@@ -377 +378 @@
       
       !                                        !* bottom temperature, salinity, velocity and depth
-      IF( kt == nit000 ) THEN
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          IF(lwp)  WRITE(numout,*) ' '
-         IF(lwp)  WRITE(numout,*) ' trabbl:bbl  : Compute bbl velocities and diffusive coefficients in ', cdtype, ' at time step ', kt
+         IF(lwp)  WRITE(numout,*) ' trabbl:bbl  : Compute bbl velocities and diffusive coefficients in ', cdtype
          IF(lwp)  WRITE(numout,*) ' '
       ENDIF
@@ -391 +392 @@
 #endif
             ik = mbkt(ji,jj)                        ! bottom T-level index
-            ztb (ji,jj) = tb(ji,jj,ik)              ! bottom before T and S
-            zsb (ji,jj) = sb(ji,jj,ik)
+            ztb (ji,jj) = tsb(ji,jj,ik,jp_tem)              ! bottom before T and S
+            zsb (ji,jj) = tsb(ji,jj,ik,jp_sal)
             zdep(ji,jj) = fsdept_0(ji,jj,ik)        ! bottom T-level reference depth
             !

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traldf.F90

@@ -62 +62 @@
       !!
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ztrdt, ztrds
-      REAL(wp), DIMENSION(jpi,jpj,jpts) :: zgtsu, zgtsv
       !!----------------------------------------------------------------------
 
@@ -70 +69 @@
       ENDIF
 
-      zgtsu(:,:,jp_tem) = gtu(:,:)        ;        zgtsu(:,:,jp_sal) = gsu(:,:)
-      zgtsv(:,:,jp_tem) = gtv(:,:)        ;        zgtsv(:,:,jp_sal) = gsv(:,:)
-
       SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend
-      CASE ( 0 )   ;   CALL tra_ldf_lap   ( kt , 'TRA', zgtsu, zgtsv, tsb, tsa, jpts        )  ! iso-level laplacian
-      CASE ( 1 )   ;   CALL tra_ldf_iso   ( kt , 'TRA', zgtsu, zgtsv, tsb, tsa, jpts, ahtb0 )  ! rotated laplacian 
-      CASE ( 2 )   ;   CALL tra_ldf_bilap ( kt , 'TRA', zgtsu, zgtsv, tsb, tsa, jpts        )  ! iso-level bilaplacian
-      CASE ( 3 )   ;   CALL tra_ldf_bilapg( kt , 'TRA',               tsb, tsa, jpts        )  ! s-coord. horizontal bilaplacian
+      CASE ( 0 )   ;   CALL tra_ldf_lap   ( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts        )  ! iso-level laplacian
+      CASE ( 1 )   ;   CALL tra_ldf_iso   ( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 )  ! rotated laplacian 
+      CASE ( 2 )   ;   CALL tra_ldf_bilap ( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts        )  ! iso-level bilaplacian
+      CASE ( 3 )   ;   CALL tra_ldf_bilapg( kt, 'TRA',             tsb, tsa, jpts        )  ! s-coord. horizontal bilaplacian
          !
       CASE ( -1 )                                     ! esopa: test all possibility with control print
-         CALL tra_ldf_lap   ( kt , 'TRA', zgtsu, zgtsv, tsb, tsa, jpts        ) 
+         CALL tra_ldf_lap   ( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts        ) 
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf0 - Ta: ', mask1=tmask,               &
          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_ldf_iso   ( kt , 'TRA', zgtsu, zgtsv, tsb, tsa, jpts, ahtb0 )  
+         CALL tra_ldf_iso   ( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 )  
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf1 - Ta: ', mask1=tmask,               &
          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_ldf_bilap ( kt , 'TRA', zgtsu, zgtsv, tsb, tsa, jpts        ) 
+         CALL tra_ldf_bilap ( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts        ) 
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf2 - Ta: ', mask1=tmask,               &
          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_ldf_bilapg( kt , 'TRA',               tsb, tsa, jpts        ) 
+         CALL tra_ldf_bilapg( kt, 'TRA',             tsb, tsa, jpts        ) 
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf3 - Ta: ', mask1=tmask,               &
          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traldf_bilap.F90

@@ -28 +28 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE diaptr          ! poleward transport diagnostics
+   USE trc_oce         ! share passive tracers/Ocean variables
 
    IMPLICIT NONE
@@ -93 +94 @@
       !!----------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_ldf_bilap : iso-level biharmonic operator'
+         IF(lwp) WRITE(numout,*) 'tra_ldf_bilap : iso-level biharmonic operator on ', cdtype
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
       ENDIF

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traldf_bilapg.F90

@@ -24 +24 @@
    USE lbclnk          ! ocean lateral boundary condition (or mpp link)
    USE diaptr          ! poleward transport diagnostics 
+   USE trc_oce         ! share passive tracers/Ocean variables
 
    IMPLICIT NONE
@@ -80 +81 @@
       !!----------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_ldf_bilapg : horizontal biharmonic operator in s-coordinate'
+         IF(lwp) WRITE(numout,*) 'tra_ldf_bilapg : horizontal biharmonic operator in s-coordinate on ', cdtype
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
       ENDIF

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traldf_iso.F90

@@ -28 +28 @@
    USE ldfslp          ! iso-neutral slopes
    USE diaptr          ! poleward transport diagnostics
+   USE trc_oce         ! share passive tracers/Ocean variables
 #if defined key_diaar5
    USE phycst          ! physical constants
@@ -116 +117 @@
       !!----------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator'
+         IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator on ', cdtype
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
       ENDIF

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traldf_lap.F90

@@ -23 +23 @@
    USE in_out_manager  ! I/O manager
    USE diaptr          ! poleward transport diagnostics
+   USE trc_oce         ! share passive tracers/Ocean variables
 
 
@@ -84 +85 @@
       !!----------------------------------------------------------------------
       
-      IF( kt == nit000 ) THEN
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_ldf_lap : iso-level laplacian diffusion'
+         IF(lwp) WRITE(numout,*) 'tra_ldf_lap : iso-level laplacian diffusion on ', cdtype
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
          e1ur(:,:) = e2u(:,:) / e1u(:,:)

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/tranpc.F90

@@ -74 +74 @@
          inpci = 0
 
-         CALL eos( tsa(:,:,:,jp_tem), tsa(:,:,:,jp_sal), rhd, zrhop )         ! Potential density
+         CALL eos( tsa, rhd, zrhop )         ! Potential density
 
          IF( l_trdtra )   THEN                    !* Save ta and sa trends

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/tranxt.F90

@@ -48 +48 @@
    PUBLIC   tra_nxt_vvl   ! to be used in trcnxt
 
-   REAL(wp), DIMENSION(jpk) ::   r2dt_t   ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler)
+   REAL(wp), DIMENSION(jpk) ::   r2dt   ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler)
 
    !! * Substitutions
@@ -120 +120 @@
  
       ! set time step size (Euler/Leapfrog)
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt_t(:) =     rdttra(:)      ! at nit000             (Euler)
-      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt_t(:) = 2.* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
+      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt(:) =     rdttra(:)      ! at nit000             (Euler)
+      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt(:) = 2.* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
       ENDIF
 
@@ -145 +145 @@
       IF( l_trdtra ) THEN              ! trend of the Asselin filter (tb filtered - tb)/dt     
          DO jk = 1, jpkm1
-            zfact = 1.e0 / r2dt_t(jk)             
+            zfact = 1.e0 / r2dt(jk)             
             ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact
             ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact
@@ -194 +194 @@
       !!----------------------------------------------------------------------
 
-      IF( kt == kit000 ) THEN
+      IF( kt == kit000 )  THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping'

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trazdf.F90

@@ -79 +79 @@
 
       SELECT CASE ( nzdf )                       ! compute lateral mixing trend and add it to the general trend
-      CASE ( 0 )    ;    CALL tra_zdf_exp( kt , 'TRA', r2dt, nn_zdfexp, tsb, tsa, jpts )  !   explicit scheme 
-      CASE ( 1 )    ;    CALL tra_zdf_imp( kt , 'TRA', r2dt,            tsb, tsa, jpts )  !   implicit scheme 
+      CASE ( 0 )    ;    CALL tra_zdf_exp( kt, 'TRA', r2dt, nn_zdfexp, tsb, tsa, jpts )  !   explicit scheme 
+      CASE ( 1 )    ;    CALL tra_zdf_imp( kt, 'TRA', r2dt,            tsb, tsa, jpts )  !   implicit scheme 
       CASE ( -1 )                                       ! esopa: test all possibility with control print
-         CALL tra_zdf_exp( kt , 'TRA', r2dt, nn_zdfexp, tsb, tsa, jpts )
+         CALL tra_zdf_exp( kt, 'TRA', r2dt, nn_zdfexp, tsb, tsa, jpts )
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' zdf0 - Ta: ', mask1=tmask,               &
          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_zdf_imp( kt , 'TRA', r2dt,            tsb, tsa, jpts ) 
+         CALL tra_zdf_imp( kt, 'TRA', r2dt,            tsb, tsa, jpts ) 
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' zdf1 - Ta: ', mask1=tmask,               &
          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trazdf_exp.F90

@@ -29 +29 @@
    USE zdfddm          ! ocean vertical physics: double diffusion
    USE in_out_manager  ! I/O manager
+   USE trc_oce         ! share passive tracers/Ocean variables
 
    IMPLICIT NONE
@@ -87 +88 @@
       !!---------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_zdf_exp : explicit vertical mixing'
+         IF(lwp) WRITE(numout,*) 'tra_zdf_exp : explicit vertical mixing on ', cdtype
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
       ENDIF

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trazdf_imp.F90

@@ -31 +31 @@
    USE domvvl          ! variable volume
    USE ldftra          ! lateral mixing type
+   USE trc_oce         ! share passive tracers/Ocean variables
 
    IMPLICIT NONE
@@ -105 +106 @@
       !!---------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          IF(lwp)WRITE(numout,*)
-         IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing'
+         IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing on ', cdtype
          IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ '
          zavi = 0.e0      ! avoid warning at compilation phase when lk_ldfslp=F

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/zpshde.F90

@@ -14 +14 @@
    !!   zps_hde      :  Horizontal DErivative of T, S and rd at the last
    !!                   ocean level (Z-coord. with Partial Steps)
-   !!   zps_hde_trc  :  Horizontal DErivative of passive tracers at the last
-   !!                   ocean level (Z-coord. with Partial Steps)
    !!----------------------------------------------------------------------
    !! * Modules used
@@ -31 +29 @@
    PUBLIC zps_hde          ! routine called by step.F90
    PUBLIC zps_hde_init     ! routine called by opa.F90
-#if defined key_top
-   PUBLIC zps_hde_trc 
-#endif
 
    !! * module variables
@@ -49 +44 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE zps_hde ( kt, ptem, psal, prd ,   &
-                            pgtu, pgsu, pgru,   &
-                            pgtv, pgsv, pgrv  )
+
+   SUBROUTINE zps_hde( kt, kjpt, pta, pgtu, pgtv,   &
+                                 prd, pgru, pgrv    )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE zps_hde  ***
@@ -91 +86 @@
       !!          di(rho) = rd~ - rd(i,j,k) or rd (i+1,j,k) - rd~
       !!
-      !! ** Action  : - pgtu, pgsu, pgru: horizontal gradient of T, S
-      !!                and rd at U-points 
-      !!              - pgtv, pgsv, pgrv: horizontal gradient of T, S
+      !! ** Action  : - pgtu, pgtv: horizontal gradient of tracer at U/V-points
+      !!              - pgru, pgrv: horizontal gradient of rd if present at U/V-points 
       !!                and rd at V-points 
       !!----------------------------------------------------------------------
       !! * Arguments
-      INTEGER, INTENT( in ) ::   kt ! ocean time-step index
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in )  :: ptem, psal, prd  ! 3D T, S and rd fields
-      REAL(wp), DIMENSION(jpi,jpj)    , INTENT( out ) :: pgtu, pgsu, pgru !  horizontal grad. of T, S and rd at u-point 
-      REAL(wp), DIMENSION(jpi,jpj)    , INTENT( out ) :: pgtv, pgsv, pgrv !  horizontal grad. of T, S and rd at v-point 
-      !! * Local declarations
-      INTEGER ::   ji , jj          ! Dummy loop indices
-      INTEGER ::   iku, ikv         ! partial step level at u- and v-points
-      REAL(wp), DIMENSION(jpi,jpj) ::   zti, ztj, zsi, zsj   ! interpolated value of T, S
-      REAL(wp), DIMENSION(jpi,jpj) ::   zri, zrj             ! interpolated value of rd
-      REAL(wp), DIMENSION(jpi,jpj) ::   zhgi, zhgj           ! depth of interpolation for eos2d
-      REAL(wp) :: ze3wu, ze3wv, zmaxu, zmaxv                 ! temporary scalars
-
-
-      ! Interpolation of T and S at the last ocean level
-# if defined key_vectopt_loop
-         jj = 1
-         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
-# else
-      DO jj = 1, jpjm1
-         DO ji = 1, jpim1
-# endif
-            ! last level
-            iku = mbatu(ji,jj)
-            ikv = mbatv(ji,jj)
-            ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
-            ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
-
-            ! i- direction
-            IF( ze3wu >= 0. ) THEN      ! case 1
-               ! interpolated values of T and S
-               zmaxu =  ze3wu / fse3w(ji+1,jj,iku)
-               zti(ji,jj) = ptem(ji+1,jj,iku) + zmaxu * ( ptem(ji+1,jj,iku-1) - ptem(ji+1,jj,iku) )
-               zsi(ji,jj) = psal(ji+1,jj,iku) + zmaxu * ( psal(ji+1,jj,iku-1) - psal(ji+1,jj,iku) )
-               ! depth of the partial step level
-               zhgi(ji,jj) = fsdept(ji,jj,iku)
-               ! gradient of T and S
-               pgtu(ji,jj) = umask(ji,jj,1) * ( zti(ji,jj) - ptem(ji,jj,iku) )
-               pgsu(ji,jj) = umask(ji,jj,1) * ( zsi(ji,jj) - psal(ji,jj,iku) )
-
-            ELSE                        ! case 2
-               ! interpolated values of T and S
-               zmaxu =  -ze3wu / fse3w(ji,jj,iku)
-               zti(ji,jj) = ptem(ji,jj,iku) + zmaxu * ( ptem(ji,jj,iku-1) - ptem(ji,jj,iku) )
-               zsi(ji,jj) = psal(ji,jj,iku) + zmaxu * ( psal(ji,jj,iku-1) - psal(ji,jj,iku) )
-               ! depth of the partial step level
-               zhgi(ji,jj) = fsdept(ji+1,jj,iku)
-               ! gradient of T and S 
-               pgtu(ji,jj) = umask(ji,jj,1) * ( ptem(ji+1,jj,iku) - zti (ji,jj) )
-               pgsu(ji,jj) = umask(ji,jj,1) * ( psal(ji+1,jj,iku) - zsi (ji,jj) )
-            ENDIF
-
-            ! j- direction
-            IF( ze3wv >= 0. ) THEN      ! case 1
-               ! interpolated values of T and S
-               zmaxv =  ze3wv / fse3w(ji,jj+1,ikv)
-               ztj(ji,jj) = ptem(ji,jj+1,ikv) + zmaxv * ( ptem(ji,jj+1,ikv-1) - ptem(ji,jj+1,ikv) )
-               zsj(ji,jj) = psal(ji,jj+1,ikv) + zmaxv * ( psal(ji,jj+1,ikv-1) - psal(ji,jj+1,ikv) )
-               ! depth of the partial step level
-               zhgj(ji,jj) = fsdept(ji,jj,ikv) 
-               ! gradient of T and S
-               pgtv(ji,jj) = vmask(ji,jj,1) * ( ztj(ji,jj) - ptem(ji,jj,ikv) )
-               pgsv(ji,jj) = vmask(ji,jj,1) * ( zsj(ji,jj) - psal(ji,jj,ikv) )
-
-            ELSE                        ! case 2
-               ! interpolated values of T and S
-               zmaxv =  -ze3wv / fse3w(ji,jj,ikv)
-               ztj(ji,jj) = ptem(ji,jj,ikv) + zmaxv * ( ptem(ji,jj,ikv-1) - ptem(ji,jj,ikv) )
-               zsj(ji,jj) = psal(ji,jj,ikv) + zmaxv * ( psal(ji,jj,ikv-1) - psal(ji,jj,ikv) ) 
-               ! depth of the partial step level
-               zhgj(ji,jj) = fsdept(ji,jj+1,ikv) 
-               ! gradient of T and S
-               pgtv(ji,jj) = vmask(ji,jj,1) * ( ptem(ji,jj+1,ikv) - ztj(ji,jj) )
-               pgsv(ji,jj) = vmask(ji,jj,1) * ( psal(ji,jj+1,ikv) - zsj(ji,jj) )
-            ENDIF
-# if ! defined key_vectopt_loop
-         END DO
-# endif
-      END DO
-
-      ! Compute interpolated rd from zti, zsi, ztj, zsj for the 2 cases at the depth of the partial
-      ! step and store it in  zri, zrj for each  case
-      CALL eos( zti, zsi, zhgi, zri )
-      CALL eos( ztj, zsj, zhgj, zrj )
-
-
-      ! Gradient of density at the last level 
-# if defined key_vectopt_loop
-         jj = 1
-         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
-# else
-      DO jj = 1, jpjm1
-         DO ji = 1, jpim1
-# endif
-            iku = mbatu(ji,jj)
-            ikv = mbatv(ji,jj)
-            ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
-            ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
-            IF( ze3wu >= 0. ) THEN    ! i-direction: case 1
-               pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji,jj) - prd(ji,jj,iku) )
-            ELSE                      ! i-direction: case 2
-               pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj) )
-            ENDIF
-            IF( ze3wv >= 0. ) THEN    ! j-direction: case 1
-               pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj) - prd(ji,jj,ikv) )  
-            ELSE                      ! j-direction: case 2
-               pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )
-            ENDIF
-# if ! defined key_vectopt_loop
-         END DO
-# endif
-      END DO
-
-      ! Lateral boundary conditions on each gradient
-      CALL lbc_lnk( pgtu , 'U', -1. )   ;   CALL lbc_lnk( pgtv , 'V', -1. )
-      CALL lbc_lnk( pgsu , 'U', -1. )   ;   CALL lbc_lnk( pgsv , 'V', -1. )
-      CALL lbc_lnk( pgru , 'U', -1. )   ;   CALL lbc_lnk( pgrv , 'V', -1. )
-
-   END SUBROUTINE zps_hde
-
-#if defined key_top
-   !!----------------------------------------------------------------------
-   !!   'key_top'                                                TOP models
-   !!----------------------------------------------------------------------
-   SUBROUTINE zps_hde_trc ( kt, kjpt, ptra, pgtru, pgtrv )
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE zps_hde_trc  ***
-      !!                    
-      !! ** Purpose :   Compute the horizontal derivative of passive tracers
-      !!      TRA at u- and v-points with a linear interpolation for z-coordinate
-      !!      with partial steps.
-      !!
-      !! ** Method  :   the same for T & S
-      !!
-      !! ** Action  : - pgtru : horizontal gradient of TRA at U-points 
-      !!              - pgtrv : horizontal gradient of TRA at V-points 
-      !!----------------------------------------------------------------------
-      !! * Arguments
-      INTEGER                              , INTENT( in )  ::  kt    ! ocean time-step index
-      INTEGER                              , INTENT( in )  ::  kjpt  ! number of tracers
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT( in )  ::  ptra  ! 4D tracers fields
-      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT( out ) ::  pgtru, pgtrv  ! horizontal grad. of TRA u- and v-points 
+      INTEGER                              , INTENT( in )           ::  kt    ! ocean time-step index
+      INTEGER                              , INTENT( in )           ::  kjpt  ! number of tracers
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT( in )           ::  pta   ! 4D active or passive tracers fields
+      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT( out)           ::  pgtu, pgtv  ! horizontal grad. of ptra u- and v-points 
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT( in ), OPTIONAL ::  prd   ! 3D rd fields
+      REAL(wp), DIMENSION(jpi,jpj         ), INTENT( out), OPTIONAL ::  pgru, pgrv  ! horizontal grad. of prd u- and v-points 
       !! * Local declarations
       INTEGER  ::   ji, jj, jn      ! Dummy loop indices
       INTEGER  ::   iku, ikv        ! partial step level at u- and v-points
-      REAL(wp) ::   ztrai, ztraj, ze3wu, ze3wv, zmaxu, zmaxv  ! temporary scalars
-      !!----------------------------------------------------------------------
-
+      REAL(wp), DIMENSION(jpi,jpj,kjpt) ::   zti, ztj     ! interpolated value of tracer
+      REAL(wp), DIMENSION(jpi,jpj)      ::   zri, zrj     ! interpolated value of rd
+      REAL(wp), DIMENSION(jpi,jpj)      ::   zhi, zhj     ! depth of interpolation for eos2d
+      REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv  ! temporary scalars
+      !!----------------------------------------------------------------------
+
+
+      ! Interpolation of tracers at the last ocean level
       DO jn = 1, kjpt
-         ! Interpolation of passive tracers at the last ocean level
 # if defined key_vectopt_loop
          jj = 1
@@ -262 +125 @@
                IF( ze3wu >= 0. ) THEN      ! case 1
                   zmaxu =  ze3wu / fse3w(ji+1,jj,iku)
-                  ! interpolated values of passive tracers
-                  ztrai = ptra(ji+1,jj,iku,jn) + zmaxu * ( ptra(ji+1,jj,iku-1,jn) - ptra(ji+1,jj,iku,jn) )
-                  ! gradient of passive tracers
-                  pgtru(ji,jj,jn) = umask(ji,jj,1) * ( ztrai - ptra(ji,jj,iku,jn) )
+                  ! interpolated values of tracers
+                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku-1,jn) - pta(ji+1,jj,iku,jn) )
+                  ! gradient of  tracers
+                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
                ELSE                        ! case 2
                   zmaxu = -ze3wu / fse3w(ji,jj,iku)
-                  ! interpolated values of passive tracers
-                  ztrai = ptra(ji,jj,iku,jn) + zmaxu * ( ptra(ji,jj,iku-1,jn) - ptra(ji,jj,iku,jn) )
-                  ! gradient of passive tracers
-                  pgtru(ji,jj,jn) = umask(ji,jj,1) * ( ptra(ji+1,jj,iku,jn) - ztrai )
+                  ! interpolated values of tracers
+                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku-1,jn) - pta(ji,jj,iku,jn) )
+                  ! gradient of tracers
+                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
                ENDIF
 
@@ -277 +140 @@
                IF( ze3wv >= 0. ) THEN      ! case 1
                   zmaxv =  ze3wv / fse3w(ji,jj+1,ikv)
-                  ! interpolated values of passive tracers
-                  ztraj = ptra(ji,jj+1,ikv,jn) + zmaxv * ( ptra(ji,jj+1,ikv-1,jn) - ptra(ji,jj+1,ikv,jn) )
-                  ! gradient of passive tracers
-                  pgtrv(ji,jj,jn) = vmask(ji,jj,1) * ( ztraj - ptra(ji,jj,ikv,jn) )
+                  ! interpolated values of tracers
+                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv-1,jn) - pta(ji,jj+1,ikv,jn) )
+                  ! gradient of tracers
+                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
                ELSE                        ! case 2
                   zmaxv =  -ze3wv / fse3w(ji,jj,ikv)
-                  ! interpolated values of passive tracers
-                  ztraj = ptra(ji,jj,ikv,jn) + zmaxv * ( ptra(ji,jj,ikv-1,jn) - ptra(ji,jj,ikv,jn) )
-                  ! gradient of passive tracers
-                  pgtrv(ji,jj,jn) = vmask(ji,jj,1) * ( ptra(ji,jj+1,ikv,jn) - ztraj )
+                  ! interpolated values of tracers
+                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv-1,jn) - pta(ji,jj,ikv,jn) )
+                  ! gradient of tracers
+                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
                ENDIF
 # if ! defined key_vectopt_loop
@@ -294 +157 @@
 
          ! Lateral boundary conditions on each gradient
-         CALL lbc_lnk( pgtru(:,:,jn) , 'U', -1. )
-         CALL lbc_lnk( pgtrv(:,:,jn) , 'V', -1. )
+         CALL lbc_lnk( pgtu(:,:,jn) , 'U', -1. )
+         CALL lbc_lnk( pgtv(:,:,jn) , 'V', -1. )
 
       END DO
 
-   END SUBROUTINE zps_hde_trc
-#endif
+      ! horizontal derivative of rd
+      IF( PRESENT( prd ) ) THEN
+         ! depth of the partial step level
+# if defined key_vectopt_loop
+         jj = 1
+         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
+# else
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+# endif
+               iku = mbatu(ji,jj)
+               ikv = mbatv(ji,jj)
+               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
+               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
+               IF( ze3wu >= 0. ) THEN   
+                  zhi(ji,jj) = fsdept(ji  ,jj,iku)
+               ELSE                   
+                  zhi(ji,jj) = fsdept(ji+1,jj,iku)
+               ENDIF
+               IF( ze3wv >= 0. ) THEN  
+                  zhj(ji,jj) = fsdept(ji,jj  ,ikv)
+               ELSE                   
+                  zhj(ji,jj) = fsdept(ji,jj+1,ikv)
+               ENDIF
+# if ! defined key_vectopt_loop
+            END DO
+# endif
+         END DO
+
+         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial
+         ! step and store it in  zri, zrj for each  case
+         CALL eos( zti, zhi, zri )
+         CALL eos( ztj, zhj, zrj )
+
+         ! Gradient of density at the last level 
+# if defined key_vectopt_loop
+         jj = 1
+         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
+# else
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+# endif
+               iku = mbatu(ji,jj)
+               ikv = mbatv(ji,jj)
+               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
+               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
+               IF( ze3wu >= 0. ) THEN    ! i-direction: case 1
+                  pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji,jj) - prd(ji,jj,iku) )
+               ELSE                      ! i-direction: case 2
+                  pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj) )
+               ENDIF
+               IF( ze3wv >= 0. ) THEN    ! j-direction: case 1
+                  pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj) - prd(ji,jj,ikv) )
+               ELSE                      ! j-direction: case 2
+                  pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )
+               ENDIF
+# if ! defined key_vectopt_loop
+            END DO
+# endif
+         END DO
+
+         ! Lateral boundary conditions on each gradient
+         CALL lbc_lnk( pgru , 'U', -1. )   ;   CALL lbc_lnk( pgrv , 'V', -1. )
+         !
+      END IF
+      !
+   END SUBROUTINE zps_hde
 
    SUBROUTINE zps_hde_init

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRD/trdicp.F90

@@ -370 +370 @@
          zkepe(:,:,:) = 0.e0
    
-         CALL eos( tn, sn, rhd, rhop )       ! now potential and in situ densities
+         CALL eos( tsn, rhd, rhop )       ! now potential and in situ densities
 
          ! Density flux at w-point

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/istate.F90

@@ -42 +42 @@
    USE dynspg_exp      ! pressure gradient schemes
    USE dynspg_ts       ! pressure gradient schemes
+   USE traswp          ! Swap arrays                      (tra_swp routine)
    
    IMPLICIT NONE
@@ -122 +123 @@
          ENDIF
 
-         CALL eos( tb, sb, rhd, rhop )        ! before potential and in situ densities
+         ta(:,:,:) = 0.       ;     sa(:,:,:) = 0.
+         CALL tra_swap
+         CALL eos( tsb, rhd, rhop )        ! before potential and in situ densities
          
          IF( ln_zps .AND. .NOT. lk_c1d )   &
-            &             CALL zps_hde( nit000, tb, sb, rhd,  &  ! Partial steps: before Horizontal DErivative
-            &                                  gtu, gsu, gru, &  ! of t, s, rd at the bottom ocean level
-            &                                  gtv, gsv, grv )
+                          CALL zps_hde( nit000, 'TRA', tsb, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
+            &                                          rhd, gru , grv , jpts )      ! of t, s, rd at the last ocean level
          
       ENDIF

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/oce.F90

@@ -39 +39 @@
    !! interpolated gradient (only used in zps case)
    !! ---------------------
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   gtu, gsu, gru   !: horizontal gradient of T, S and rd at bottom u-point
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   gtv, gsv, grv   !: horizontal gradient of T, S and rd at bottom v-point 
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpts) ::   gtsu, gtsv !: horizontal gradient of T, S bottom u-point
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)      ::   gru , grv   !: horizontal gradient of rd at bottom u-point
 
    !! free surface                       !  before  !  now     !  after   !

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/opa.F90

@@ -37 +37 @@
    USE step_oce        ! Time stepping module definition
    USE sbc_oce         ! surface boundary condition: ocean
-   USE trdmod_oce      ! ocean variables trends
    USE domcfg          ! domain configuration               (dom_cfg routine)
    USE mppini          ! shared/distributed memory setting (mpp_init routine)
@@ -48 +47 @@
    USE zdfini
    USE phycst          ! physical constant                  (par_cst routine)
-   USE trdmod          ! momentum/tracers trends       (trd_mod_init routine)
    USE step            ! OPA time-stepping                  (stp     routine)
 #if defined key_oasis3
@@ -160 +158 @@
       !!
       !!----------------------------------------------------------------------
-#if defined key_oasis3 || defined key_oasis4 || defined key_iomput
       INTEGER :: ilocal_comm
-#endif
       CHARACTER(len=80),dimension(10) ::   cltxt = ''
       INTEGER                         ::   ji   ! local loop indices
@@ -192 +188 @@
       narea = mynode( cltxt, ilocal_comm )    ! Nodes selection (control print return in cltxt)
 # else
+      ilocal_comm = 0
       narea = mynode( cltxt )                 ! Nodes selection (control print return in cltxt)
 # endif
@@ -220 +217 @@
       !                             !--------------------------------!
 
-      CALL opa_flg                          ! Control prints & Benchmark
-
-                                            ! Domain decomposition
+      CALL opa_flg                           ! Control prints & Benchmark
+
+      !                                      ! Domain decomposition
       IF( jpni*jpnj == jpnij ) THEN   ;   CALL mpp_init      ! standard cutting out
       ELSE                            ;   CALL mpp_init2     ! eliminate land processors
@@ -231 +228 @@
      
 
+      !                                      ! General initialization
                             CALL     phy_cst    ! Physical constants
                             CALL     eos_init   ! Equation of state
@@ -245 +243 @@
 
       !                                     ! Ocean physics
-                            CALL     sbc_init   ! Read namsbc namelist : surface module (needed for iom_init)
-
-      !                                     ! Vertical physics
-                            CALL     zdf_init   ! namelist read
-                            CALL zdf_bfr_init   ! bottom friction
-      IF( lk_zdfric     )   CALL zdf_ric_init   ! Richardson number dependent Kz
-      IF( lk_zdftke_old )   CALL zdf_tke_init   ! TKE closure scheme for Kz (old scheme)
-      IF( lk_zdftke     )   CALL     tke_init   ! TKE closure scheme for Kz
-      IF( lk_zdfkpp     )   CALL zdf_kpp_init   ! KPP closure scheme for Kz
-      IF( lk_zdftmx     )   CALL zdf_tmx_init   ! tidal vertical mixing
+                            CALL     sbc_init   !  Forcings : surface module 
+
+      !                                         ! Vertical physics
+                            CALL     zdf_init      ! namelist read
+                            CALL zdf_bfr_init      ! bottom friction
+      IF( lk_zdfric     )   CALL zdf_ric_init      ! Richardson number dependent Kz
+      IF( lk_zdftke_old )   CALL zdf_tke_init      ! TKE closure scheme for Kz (old scheme)
+      IF( lk_zdftke     )   CALL     tke_init      ! TKE closure scheme for Kz
+      IF( lk_zdfkpp     )   CALL zdf_kpp_init      ! KPP closure scheme for Kz
+      IF( lk_zdftmx     )   CALL zdf_tmx_init      ! tidal vertical mixing
       IF( lk_zdfddm .AND. .NOT. lk_zdfkpp )   & 
-         &                  CALL zdf_ddm_init   ! double diffusive mixing
-      !                                     ! Lateral physics
-                            CALL ldf_tra_init   ! Lateral ocean tracer physics
-                            CALL ldf_dyn_init   ! Lateral ocean momentum physics
-      IF( lk_ldfslp )       CALL ldf_slp_init   ! slope of lateral mixing
+         &                  CALL zdf_ddm_init      ! double diffusive mixing
+      !                                         ! Lateral physics
+                            CALL ldf_tra_init      ! Lateral ocean tracer physics
+                            CALL ldf_dyn_init      ! Lateral ocean momentum physics
+      IF( lk_ldfslp )       CALL ldf_slp_init      ! slope of lateral mixing
+
       !                                     ! Active tracers
                             CALL tra_qsr_init   ! penetrative solar radiation qsr
@@ -280 +279 @@
 #if defined key_top
       !                                     ! Passive tracers
-                            CALL     trc_ini
-#endif
-
+                            CALL     trc_init
+#endif
       !                                     ! diagnostics
                             CALL     iom_init   ! iom_put initialization

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/step.F90

@@ -105 +105 @@
       ! Ocean physics update                (ua, va, ta, sa used as workspace)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL bn2( tb, sb, rn2b )     ! before Brunt-Vaisala frequency
-                         CALL bn2( tn, sn, rn2  )     ! now    Brunt-Vaisala frequency
+                         CALL bn2( tsb, rn2b )        ! before Brunt-Vaisala frequency
+                         CALL bn2( tsn, rn2  )        ! now    Brunt-Vaisala frequency
       !
       !  VERTICAL PHYSICS   
@@ -137 +137 @@
       !
       IF( lk_ldfslp ) THEN                            ! slope of lateral mixing
-                         CALL eos( tb, sb, rhd )                ! before in situ density
-         IF( ln_zps )    CALL zps_hde( kstp, tb, sb, rhd,  &    ! Partial steps: before horizontal gradient
-            &                              gtu, gsu, gru,  &    ! of t, s, rd at the last ocean level
-            &                              gtv, gsv, grv )
+                         CALL eos( tsb, rhd )                ! before in situ density
+         IF( ln_zps )    CALL zps_hde( kstp, jpts, tsb, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
+            &                                      rhd, gru , grv  )      ! of t, s, rd at the last ocean level
                          CALL ldf_slp( kstp, rhd, rn2b )        ! before slope of the lateral mixing
       ENDIF
@@ -169 +168 @@
       ! Active tracers                              (ua, va used as workspace)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                             ta(:,:,:) = 0.e0               ! set tracer trends to zero
-                             sa(:,:,:) = 0.e0
-
-                             CALL tra_swap
+                             tsa(:,:,:,:) = 0.e0               ! set tracer trends to zero
+
                              CALL tra_sbc    ( kstp )       ! surface boundary condition
       IF( ln_traqsr      )   CALL tra_qsr    ( kstp )       ! penetrative solar radiation qsr
@@ -192 +189 @@
          IF( ln_zdfnpc   )   CALL tra_npc    ( kstp )            ! update after fields by non-penetrative convection
                              CALL tra_nxt    ( kstp )            ! tracer fields at next time step
-                             CALL tra_unswap
-                             CALL eos( ta, sa, rhd, rhop )       ! Time-filtered in situ density for hpg computation
-         IF( ln_zps      )   CALL zps_hde( kstp, ta, sa, rhd,   &   ! Partial steps: time filtered hor. derivative
-            &                                   gtu, gsu, gru,  &   ! of t, s, rd at the bottom ocean level
-            &                                   gtv, gsv, grv )
+                             CALL eos( tsa, rhd, rhop )       ! Time-filtered in situ density for hpg computation
+         IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsa, gtsu, gtsv,  &    ! Partial steps: time filtered hor. derivative
+            &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level
          
       ELSE                                                  ! centered hpg  (eos then time stepping)
-                             CALL tra_unswap
-                             CALL eos( tn, sn, rhd, rhop )       ! now in situ density for hpg computation
-         IF( ln_zps      )   CALL zps_hde( kstp, tn, sn, rhd,   &   ! Partial steps: now horizontal derivative
-            &                                   gtu, gsu, gru,  &   ! of t, s, rd at the bottom ocean level
-            &                                   gtv, gsv, grv )
-                             CALL tra_swap
+                             CALL eos( tsn, rhd, rhop )       ! now in situ density for hpg computation
+         IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: now hor. derivative
+            &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level
          IF( ln_zdfnpc   )   CALL tra_npc    ( kstp )       ! update after fields by non-penetrative convection
                              CALL tra_nxt    ( kstp )       ! tracer fields at next time step
+      ENDIF 
                              CALL tra_unswap
-      ENDIF 
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/step_oce.F90

@@ -81 +81 @@
    USE trdmld_rst      ! restart for mixed-layer trends
    USE trdmod_oce      ! ocean momentum/tracers trends
+   USE trdmod          ! momentum/tracers trends   
    USE trdvor          ! vorticity budget                 (trd_vor routine)
    USE diagap          ! hor. mean model-data gap         (dia_gap routine)

branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/trc_oce.F90

@@ -20 +20 @@
    PUBLIC   trc_oce_rgb_read   ! routine called by traqsr.F90
    PUBLIC   trc_oce_ext_lev    ! function called by traqsr.F90 at least
-   
+ 
+   INTEGER , PUBLIC  ::  nn_dttrc      !: frequency of step on passive tracers
+   INTEGER , PUBLIC  ::  nittrc000     !: first time step of passive tracers model
+  
    REAL(wp), PUBLIC , DIMENSION(jpi,jpj,jpk) ::   etot3   !: light absortion coefficient
 

branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/C14b/trcsms_c14b.F90

@@ -234 +234 @@
             !  Computation of solubility  
             IF (tmask(ji,jj,1) >  0.) THEN
-               ztp  = ( tn(ji,jj,1) + 273.16 ) * 0.01
+               ztp  = ( tsn(ji,jj,1,jp_tem) + 273.16 ) * 0.01
                zsk  = 0.023517 + ztp * ( -0.023656 + 0.0047036 * ztp )
-               zsol = EXP( -60.2409 + 93.4517 / ztp  + 23.3585 * LOG( ztp ) + zsk * sn(ji,jj,1) )
+               zsol = EXP( -60.2409 + 93.4517 / ztp  + 23.3585 * LOG( ztp ) + zsk * tsn(ji,jj,1,jp_sal) )
                ! convert solubilities [mol/(l * atm)] -> [mol/(m^3 * ppm)]
                zsol = zsol * 1.0e-03
@@ -247 +247 @@
 
             ! Computes the Schmidt number of CO2 in seawater
-            zt   = tn(ji,jj,1)
+            zt   = tsn(ji,jj,1,jp_tem)
             zsch = 2073.1 + zt * ( -125.62 + zt * (3.6276 - 0.043219 * zt ) )
 

branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/CFC/trcsms_cfc.F90

@@ -129 +129 @@
                ! coefficient for solubility for CFC-11/12 in  mol/l/atm
                IF( tmask(ji,jj,1) .GE. 0.5 ) THEN
-                  ztap  = ( tn(ji,jj,1) + 273.16 ) * 0.01
+                  ztap  = ( tsn(ji,jj,1,jp_tem) + 273.16 ) * 0.01
                   zdtap = sob(1,jl) + ztap * ( sob(2,jl) + ztap * sob(3,jl) ) 
                   zsol  =  EXP( soa(1,jl) + soa(2,jl) / ztap + soa(3,jl) * LOG( ztap )   &
-                     &                    + soa(4,jl) * ztap * ztap + sn(ji,jj,1) * zdtap ) 
+                     &                    + soa(4,jl) * ztap * ztap + tsn(ji,jj,1,jp_sal) * zdtap ) 
                ELSE
                   zsol  = 0.e0
@@ -143 +143 @@
                ! Computation of speed transfert
                !    Schmidt number
-               zt1  = tn(ji,jj,1)
+               zt1  = tsn(ji,jj,1,jp_tem)
                zt2  = zt1 * zt1 
                zt3  = zt1 * zt2

branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/LOBSTER/trcbio.F90

@@ -81 +81 @@
       !!---------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
+      IF( kt == nittrc000 ) THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) ' trc_bio: LOBSTER bio-model'

branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/LOBSTER/trcexp.F90

@@ -60 +60 @@
       !!---------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
+      IF( kt == nittrc000 ) THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) ' trc_exp: LOBSTER export'

branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/LOBSTER/trcopt.F90

@@ -65 +65 @@
       !!---------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
+      IF( kt == nittrc000 ) THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) ' trc_opt : LOBSTER optic-model'

branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/LOBSTER/trcsed.F90

@@ -67 +67 @@
       !!---------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
+      IF( kt == nittrc000 ) THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) ' trc_sed: LOBSTER sedimentation'

branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zche.F90

@@ -181 +181 @@
 
             !                             ! SET ABSOLUTE TEMPERATURE
-            ztkel = tn(ji,jj,1) + 273.16
+            ztkel = tsn(ji,jj,1,jp_tem) + 273.16
             zqtt  = ztkel * 0.01
             zqtt2 = zqtt * zqtt
-            zsal  = sn(ji,jj,1) + (1.- tmask(ji,jj,1) ) * 35.
+            zsal  = tsn(ji,jj,1,jp_sal) + (1.- tmask(ji,jj,1) ) * 35.
             zlqtt = LOG( zqtt )
 
@@ -214 +214 @@
 
                ! SET ABSOLUTE TEMPERATURE
-               ztkel   = tn(ji,jj,jk) + 273.16
+               ztkel   = tsn(ji,jj,jk,jp_tem) + 273.16
                zqtt    = ztkel * 0.01
-               zsal    = sn(ji,jj,jk) + ( 1.-tmask(ji,jj,jk) ) * 35.
+               zsal    = tsn(ji,jj,jk,jp_sal) + ( 1.-tmask(ji,jj,jk) ) * 35.
                zsqrt  = SQRT( zsal )
                zsal15  = zsqrt * zsal
@@ -224 +224 @@
                zis2   = zis * zis
                zisqrt = SQRT( zis )
-               ztc     = tn(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 20.
+               ztc     = tsn(ji,jj,jk,jp_tem) + ( 1.- tmask(ji,jj,jk) ) * 20.
 
                ! CHLORINITY (WOOSTER ET AL., 1969)

branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zflx.F90

@@ -128 +128 @@
 !CDIR NOVERRCHK
          DO ji = 1, jpi
-            ztc  = MIN( 35., tn(ji,jj,1) )
+            ztc  = MIN( 35., tsn(ji,jj,1,jp_tem) )
             ztc2 = ztc * ztc
             ztc3 = ztc * ztc2 

branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zint.F90

@@ -55 +55 @@
       ! -------------------------------------------
 
-      tgfunc (:,:,:) = EXP( 0.063913 * tn(:,:,:) )
-      tgfunc2(:,:,:) = EXP( 0.07608  * tn(:,:,:) )
+      tgfunc (:,:,:) = EXP( 0.063913 * tsn(:,:,:,jp_tem) )
+      tgfunc2(:,:,:) = EXP( 0.07608  * tsn(:,:,:,jp_tem) )
 
       ! Computation of the silicon dependant half saturation

branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zlim.F90

@@ -161 +161 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               ztemp = MAX( 0., tn(ji,jj,jk) )
+               ztemp = MAX( 0., tsn(ji,jj,jk,jp_tem) )
                xfracal(ji,jj,jk) = caco3r * xlimphy(ji,jj,jk)   &
                   &                       * MAX( 0.0001, ztemp / ( 2.+ ztemp ) )   &

branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zprod.F90

@@ -147 +147 @@
                ! Computation of the P-I slope for nanos and diatoms
                IF( etot(ji,jj,jk) > 1.E-3 ) THEN
-                   ztn    = MAX( 0., tn(ji,jj,jk) - 15. )
+                   ztn    = MAX( 0., tsn(ji,jj,jk,jp_tem) - 15. )
                    zadap  = 0.+ 1.* ztn / ( 2.+ ztn )
                    zadap2 = 0.e0

branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zrem.F90

@@ -287 +287 @@
                zsatur  = ( sio3eq(ji,jj,jk) - trn(ji,jj,jk,jpsil) ) / ( sio3eq(ji,jj,jk) + rtrn )
                zsatur  = MAX( rtrn, zsatur )
-               zsatur2 = zsatur * ( 1. + tn(ji,jj,jk) / 400.)**4
-               znusil  = 0.225  * ( 1. + tn(ji,jj,jk) / 15.) * zsatur + 0.775 * zsatur2**9
+               zsatur2 = zsatur * ( 1. + tsn(ji,jj,jk,jp_tem) / 400.)**4
+               znusil  = 0.225  * ( 1. + tsn(ji,jj,jk,jp_tem) / 15.) * zsatur + 0.775 * zsatur2**9
 #    if defined key_degrad
                zsiremin = xsirem * xstep * znusil * facvol(ji,jj,jk)

branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zsed.F90

@@ -383 +383 @@
       imois2 = nmonth
 
-      ! 1. first call kt=nit000
+      ! 1. first call kt=nittrc000
       ! -----------------------
 
-      IF( kt == nit000 ) THEN
+      IF( kt == nittrc000 ) THEN
          ! initializations
          nflx1  = 0
@@ -402 +402 @@
       ! ----------------
 
-      IF( kt == nit000 .OR. imois /= nflx1 ) THEN
+      IF( kt == nittrc000 .OR. imois /= nflx1 ) THEN
 
          ! Calendar computation
@@ -423 +423 @@
          CALL iom_get ( numdust, jpdom_data, 'dust', dustmo(:,:,1), nflx1 )
          CALL iom_get ( numdust, jpdom_data, 'dust', dustmo(:,:,2), nflx2 )
-
-         IF(lwp .AND. nitend-nit000 <= 100 ) THEN
-            WRITE(numout,*)
-            WRITE(numout,*) ' read clio flx ok'
-            WRITE(numout,*)
-               WRITE(numout,*)
-               WRITE(numout,*) 'Clio month: ',nflx1,'  field: dust'
-               CALL prihre( dustmo(:,:,1),jpi,jpj,1,jpi,20,1,jpj,10,1e9,numout )
-         ENDIF
 
       ENDIF

branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/SED/par_sed.F90

@@ -16 +16 @@
       jpjm1    =>   jpjm1 ,  & !: jpj - 1
       jpij     =>   jpij       !: jpi x jpj
+      jp_tem   =>   jp_tem     !: indice of temperature
+      jp_sal   =>   jp_sal     !: indice of salintity
 
 #if ! defined key_sed_off

branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/SED/sed.F90

@@ -34 +34 @@
 
    USE oce , ONLY :            &
-      tn      =>    tn    ,  & !: pot. temperature (celsius)
-      sn      =>    sn         !: salinity (psu)
+      tsn      =>    tsn       & !: pot. temperature (celsius) and salinity (psu)
 
    USE trc, ONLY :  &
@@ -215 +214 @@
    INTEGER, PUBLIC :: &
      numsed = 27
-   
 #else
    !!======================================================================

branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/SED/seddta.F90

@@ -131 +131 @@
                   trc_data(ji,jj,9 ) = sinking2(ji,jj,ikt)
                   trc_data(ji,jj,10) = sinkcal (ji,jj,ikt)
-                  trc_data(ji,jj,11) = tn      (ji,jj,ikt)
-                  trc_data(ji,jj,12) = sn      (ji,jj,ikt)
+                  trc_data(ji,jj,11) = tsn     (ji,jj,ikt,jp_tem)
+                  trc_data(ji,jj,12) = tsn     (ji,jj,ikt,jp_sal)
 #   else
                   trc_data(ji,jj,7 ) = sinksil (ji,jj,ikt)
                   trc_data(ji,jj,8 ) = sinking (ji,jj,ikt)
                   trc_data(ji,jj,9 ) = sinkcal (ji,jj,ikt)
-                  trc_data(ji,jj,10) = tn      (ji,jj,ikt)
-                  trc_data(ji,jj,11) = sn      (ji,jj,ikt)       
+                  trc_data(ji,jj,10) = tsn     (ji,jj,ikt,jp_tem)
+                  trc_data(ji,jj,11) = tsn     (ji,jj,ikt,jp_sal)       
 #   endif
                ENDIF

branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trcadv.F90

@@ -33 +33 @@
 
    PUBLIC   trc_adv    ! routine called by step module
+
    INTEGER ::   nadv   ! choice of the type of advection scheme
+   REAL(wp), DIMENSION(jpk) ::   r2dt  ! vertical profile time-step, = 2 rdttra
+      !                                ! except at nit000 (=rdttra) if neuler=0
 
    !! * Substitutions
@@ -64 +67 @@
       IF( kt == nittrc000 )   CALL trc_adv_ctl          ! initialisation & control of options
 
+      IF( neuler == 0 .AND. kt == nittrc000 ) THEN     ! at nit000
+         r2dt(:) =  rdttra(:) * FLOAT(nn_dttrc)          ! = rdtra (restarting with Euler time stepping)
+      ELSEIF( kt <= nittrc000 + nn_dttrc ) THEN          ! at nit000 or nit000+1
+         r2dt(:) = 2. * rdttra(:) * FLOAT(nn_dttrc)      ! = 2 rdttra (leapfrog)
+      ENDIF
+
       !                                                   ! effective transport
       DO jk = 1, jpkm1
@@ -78 +87 @@
       !
       SELECT CASE ( nadv )                            !==  compute advection trend and add it to general trend  ==!
-      CASE ( 1 )   ;    CALL tra_adv_cen2  ( kt, 'TRC', zun, zvn, zwn, trb, trn, tra, jptra )   !  2nd order centered
-      CASE ( 2 )   ;    CALL tra_adv_tvd   ( kt, 'TRC', zun, zvn, zwn, trb, trn, tra, jptra )   !  TVD 
-      CASE ( 3 )   ;    CALL tra_adv_muscl ( kt, 'TRC', zun, zvn, zwn, trb,      tra, jptra )   !  MUSCL 
-      CASE ( 4 )   ;    CALL tra_adv_muscl2( kt, 'TRC', zun, zvn, zwn, trb, trn, tra, jptra )   !  MUSCL2 
-      CASE ( 5 )   ;    CALL tra_adv_ubs   ( kt, 'TRC', zun, zvn, zwn, trb, trn, tra, jptra )   !  UBS 
-      CASE ( 6 )   ;    CALL tra_adv_qck   ( kt, 'TRC', zun, zvn, zwn, trb, trn, tra, jptra )   !  QUICKEST 
+      CASE ( 1 )   ;    CALL tra_adv_cen2  ( kt, 'TRC',       zun, zvn, zwn, trb, trn, tra, jptra )   !  2nd order centered
+      CASE ( 2 )   ;    CALL tra_adv_tvd   ( kt, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  TVD 
+      CASE ( 3 )   ;    CALL tra_adv_muscl ( kt, 'TRC', r2dt, zun, zvn, zwn, trb,      tra, jptra )   !  MUSCL 
+      CASE ( 4 )   ;    CALL tra_adv_muscl2( kt, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  MUSCL2 
+      CASE ( 5 )   ;    CALL tra_adv_ubs   ( kt, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  UBS 
+      CASE ( 6 )   ;    CALL tra_adv_qck   ( kt, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  QUICKEST 
       !
       CASE (-1 )                                      !==  esopa: test all possibility with control print  ==!
-         CALL tra_adv_cen2  ( kt, 'TRC', zun, zvn, zwn, trb, trn, tra, jptra )          
+         CALL tra_adv_cen2  ( kt, 'TRC',       zun, zvn, zwn, trb, trn, tra, jptra )          
          WRITE(charout, FMT="('adv1')")  ; CALL prt_ctl_trc_info(charout)
                                            CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
-         CALL tra_adv_tvd   ( kt, 'TRC', zun, zvn, zwn, trb, trn, tra, jptra )          
+         CALL tra_adv_tvd   ( kt, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )          
          WRITE(charout, FMT="('adv2')")  ; CALL prt_ctl_trc_info(charout)
                                            CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
-         CALL tra_adv_muscl ( kt, 'TRC', zun, zvn, zwn, trb,      tra, jptra )          
+         CALL tra_adv_muscl ( kt, 'TRC', r2dt, zun, zvn, zwn, trb,      tra, jptra )          
          WRITE(charout, FMT="('adv3')")  ; CALL prt_ctl_trc_info(charout)
                                            CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
-         CALL tra_adv_muscl2( kt, 'TRC', zun, zvn, zwn, trb, trn, tra, jptra )          
+         CALL tra_adv_muscl2( kt, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )          
          WRITE(charout, FMT="('adv4')")  ; CALL prt_ctl_trc_info(charout)
                                            CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
-         CALL tra_adv_ubs   ( kt, 'TRC', zun, zvn, zwn, trb, trn, tra, jptra )          
+         CALL tra_adv_ubs   ( kt, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )          
          WRITE(charout, FMT="('adv5')")  ; CALL prt_ctl_trc_info(charout)
                                            CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
-         CALL tra_adv_qck   ( kt, 'TRC', zun, zvn, zwn, trb, trn, tra, jptra )          
+         CALL tra_adv_qck   ( kt, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )          
          WRITE(charout, FMT="('adv6')")  ; CALL prt_ctl_trc_info(charout)
                                            CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')

branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trcldf.F90

@@ -70 +70 @@
 
       SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend
-      CASE ( 0 )   ;   CALL tra_ldf_lap   ( kt , 'TRC', gtru, gtrv, trb, tra, jptra            )  ! iso-level laplacian
-      CASE ( 1 )   ;   CALL tra_ldf_iso   ( kt , 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 )  ! rotated laplacian 
-      CASE ( 2 )   ;   CALL tra_ldf_bilap ( kt , 'TRC', gtru, gtrv, trb, tra, jptra            )  ! iso-level bilaplacian
-      CASE ( 3 )   ;   CALL tra_ldf_bilapg( kt , 'TRC',             trb, tra, jptra            )  ! s-coord. horizontal bilaplacian
+      CASE ( 0 )   ;   CALL tra_ldf_lap   ( kt, 'TRC', gtru, gtrv, trb, tra, jptra            )  ! iso-level laplacian
+      CASE ( 1 )   ;   CALL tra_ldf_iso   ( kt, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 )  ! rotated laplacian 
+      CASE ( 2 )   ;   CALL tra_ldf_bilap ( kt, 'TRC', gtru, gtrv, trb, tra, jptra            )  ! iso-level bilaplacian
+      CASE ( 3 )   ;   CALL tra_ldf_bilapg( kt, 'TRC',             trb, tra, jptra            )  ! s-coord. horizontal bilaplacian
          !
       CASE ( -1 )                                     ! esopa: test all possibility with control print
-         CALL tra_ldf_lap   ( kt , 'TRC', gtru, gtrv, trb, tra, jptra            )
+         CALL tra_ldf_lap   ( kt, 'TRC', gtru, gtrv, trb, tra, jptra            )
          WRITE(charout, FMT="('ldf0 ')") ;  CALL prt_ctl_trc_info(charout)
                                             CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
-         CALL tra_ldf_iso   ( kt , 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 )
+         CALL tra_ldf_iso   ( kt, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 )
          WRITE(charout, FMT="('ldf1 ')") ;  CALL prt_ctl_trc_info(charout)
                                             CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
-         CALL tra_ldf_bilap ( kt , 'TRC', gtru, gtrv, trb, tra, jptra            )
+         CALL tra_ldf_bilap ( kt, 'TRC', gtru, gtrv, trb, tra, jptra            )
          WRITE(charout, FMT="('ldf2 ')") ;  CALL prt_ctl_trc_info(charout)
                                             CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
-         CALL tra_ldf_bilapg( kt , 'TRC',             trb, tra, jptra            )
+         CALL tra_ldf_bilapg( kt, 'TRC',             trb, tra, jptra            )
          WRITE(charout, FMT="('ldf3 ')") ;  CALL prt_ctl_trc_info(charout)
                                             CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )

branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trcnxt.F90

@@ -44 +44 @@
    PUBLIC trc_nxt          ! routine called by step.F90
 
-  REAL(wp), DIMENSION(jpk) ::   r2dt_t
+  REAL(wp), DIMENSION(jpk) ::   r2dt
    !!----------------------------------------------------------------------
    !!   TOP 1.0 , LOCEAN-IPSL (2005) 
@@ -109 +109 @@
 
       ! set time step size (Euler/Leapfrog)
-      IF( neuler == 0 .AND. kt ==  nittrc000) THEN  ;  r2dt_t(:) =     rdttra(:) * FLOAT( nn_dttrc )  ! at nit000             (Euler)
-      ELSEIF( kt <= nittrc000 + 1 )           THEN  ;  r2dt_t(:) = 2.* rdttra(:) * FLOAT( nn_dttrc )  ! at nit000 or nit000+1 (Leapfrog)
+      IF( neuler == 0 .AND. kt ==  nittrc000) THEN  ;  r2dt(:) =     rdttra(:) * FLOAT( nn_dttrc )  ! at nit000             (Euler)
+      ELSEIF( kt <= nittrc000 + 1 )           THEN  ;  r2dt(:) = 2.* rdttra(:) * FLOAT( nn_dttrc )  ! at nit000 or nit000+1 (Leapfrog)
       ENDIF
 
@@ -133 +133 @@
          DO jn = 1, jptra
             DO jk = 1, jpkm1
-               zfact = 1.e0 / r2dt_t(jk)  
+               zfact = 1.e0 / r2dt(jk)  
                ztrdt(:,:,jk,jn) = ( trb(:,:,jk,jn) - ztrdt(:,:,jk,jn) ) * zfact 
                CALL trd_tra( kt, 'TRC', jn, jptra_trd_atf, ztrdt )

branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trctrp.F90

@@ -28 +28 @@
    USE trcrad          ! positivity                          (trc_rad routine)
    USE trcsbc          ! surface boundary condition          (trc_sbc routine)
-   USE zpshde          ! partial step: hor. derivative   (zps_hde_trc routine)
+   USE zpshde          ! partial step: hor. derivative       (zps_hde routine)
 
 #if defined key_agrif
@@ -74 +74 @@
 #endif
                                 CALL trc_zdf( kstp )            ! vertical mixing and after tracer fields
-         IF( ln_zps    )        CALL zps_hde_trc( kstp, jptra, trn, gtru, gtrv )  ! Partial steps: now horizontal gradient of passive
+         IF( ln_zps    )        CALL zps_hde( kstp, jptra, trn, gtru, gtrv )  ! Partial steps: now horizontal gradient of passive
                                                                 ! tracers at the bottom ocean level
                                 CALL trc_nxt( kstp )            ! tracer fields at next time step     

branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trczdf.F90

@@ -80 +80 @@
       SELECT CASE ( nzdf )                       ! compute lateral mixing trend and add it to the general trend
       CASE ( -1 )                                       ! esopa: test all possibility with control print
-         CALL tra_zdf_exp( kt , 'TRC', r2dt, nn_trczdf_exp, trb, tra, jptra ) 
+         CALL tra_zdf_exp( kt, 'TRC', r2dt, nn_trczdf_exp, trb, tra, jptra ) 
          WRITE(charout, FMT="('zdf1 ')") ;  CALL prt_ctl_trc_info(charout)
                                             CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
-         CALL tra_zdf_imp( kt , 'TRC', r2dt,                trb, tra, jptra ) 
+         CALL tra_zdf_imp( kt, 'TRC', r2dt,                trb, tra, jptra ) 
          WRITE(charout, FMT="('zdf2 ')") ;  CALL prt_ctl_trc_info(charout)
                                             CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
-      CASE ( 0 ) ;  CALL tra_zdf_exp( kt , 'TRC', r2dt, nn_trczdf_exp, trb, tra, jptra )    !   explicit scheme 
-      CASE ( 1 ) ;  CALL tra_zdf_imp( kt , 'TRC', r2dt,                trb, tra, jptra )    !   implicit scheme          
+      CASE ( 0 ) ;  CALL tra_zdf_exp( kt, 'TRC', r2dt, nn_trczdf_exp, trb, tra, jptra )    !   explicit scheme 
+      CASE ( 1 ) ;  CALL tra_zdf_imp( kt, 'TRC', r2dt,                trb, tra, jptra )    !   implicit scheme          
 
       END SELECT

branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/oce_trc.F90

@@ -173 +173 @@
    USE oce , ONLY :   l_traldf_rot => l_traldf_rot  !: rotated laplacian operator for lateral diffusion
 #if defined key_offline
-   USE oce , ONLY :   gtu     =>    gtu     !: t-, s- and rd horizontal gradient at u- and
-   USE oce , ONLY :   gsu     =>    gsu     !: v-points at bottom ocean level
+   USE oce , ONLY :   gtsu    =>    gtsu    !: t-, s- and rd horizontal gradient at u- and
+   USE oce , ONLY :   gtsv    =>    gtsv    !:
    USE oce , ONLY :   gru     =>    gru     !:
-   USE oce , ONLY :   gtv     =>    gtv     !:
-   USE oce , ONLY :   gsv     =>    gsv     !:
    USE oce , ONLY :   grv     =>    grv     !: 
 # if defined key_degrad

branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/trc.F90

@@ -55 +55 @@
    !! passive tracers restart (input and output)
    !! ------------------------------------------  
-   INTEGER , PUBLIC  ::  nn_dttrc      !: frequency of step on passive tracers
-   INTEGER , PUBLIC  ::  nittrc000     !: first time step of passive tracers model
    LOGICAL , PUBLIC  ::  ln_rsttr      !: boolean term for restart i/o for passive tracers (namelist)
    LOGICAL , PUBLIC  ::  lrst_trc      !: logical to control the trc restart write

branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/trcini.F90

@@ -29 +29 @@
    USE daymod
 #endif
-   USE zpshde          ! partial step: hor. derivative   (zps_hde_trc routine)
+   USE zpshde          ! partial step: hor. derivative   (zps_hde routine)
    USE in_out_manager  ! I/O manager
    USE prtctl_trc      ! Print control passive tracers (prt_ctl_trc_init routine)
@@ -49 +49 @@
 CONTAINS
    
-   SUBROUTINE trc_ini
+   SUBROUTINE trc_init
       !!---------------------------------------------------------------------
-      !!                     ***  ROUTINE trc_ini  ***
+      !!                     ***  ROUTINE trc_init  ***
       !!
       !! ** Purpose :   Initialization of the passive tracer fields 
@@ -67 +67 @@
 
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'trc_ini : initial set up of the passive tracers'
+      IF(lwp) WRITE(numout,*) 'trc_init : initial set up of the passive tracers'
       IF(lwp) WRITE(numout,*) '~~~~~~~'
 
@@ -138 +138 @@
       
       IF( ln_zps .AND. .NOT. lk_trc_c1d )   &              ! Partial steps: before horizontal gradient of passive
-      &                     CALL zps_hde_trc( nittrc000, jptra, trb, gtru, gtrv )       ! tracers at the bottom ocean level
+      &                     CALL zps_hde( nittrc000, jptra, trb, gtru, gtrv )       ! tracers at the bottom ocean level
 
 
@@ -170 +170 @@
       ENDIF
 
-   END SUBROUTINE trc_ini
+   END SUBROUTINE trc_init
 
 #else
@@ -177 +177 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE trc_ini                      ! Dummy routine   
-   END SUBROUTINE trc_ini
+   SUBROUTINE trc_init                      ! Dummy routine   
+   END SUBROUTINE trc_init
 #endif