Changeset 2148

Timestamp:

2010-10-04T15:53:42+02:00 (14 years ago)

Author:

cetlod

Message:

merge LOCEAN 2010 developments branches

Location:

branches/DEV_r2106_LOCEAN2010/NEMO

Files:

• OPA_SRC/DIA/diahsb.F90 (added)
• OPA_SRC/DOM/dom_oce.F90 (modified) (1 diff)
• OPA_SRC/DOM/domvvl.F90 (modified) (3 diffs)
• OPA_SRC/DOM/domzgr_substitute.h90 (modified) (4 diffs)
• OPA_SRC/DYN/divcur.F90 (modified) (3 diffs)
• OPA_SRC/DYN/dynnxt.F90 (modified) (6 diffs)
• OPA_SRC/DYN/dynzdf_exp.F90 (modified) (8 diffs)
• OPA_SRC/DYN/dynzdf_imp.F90 (modified) (17 diffs)
• OPA_SRC/DYN/sshwzv.F90 (modified) (15 diffs)
• OPA_SRC/IOM/restart.F90 (modified) (7 diffs)
• OPA_SRC/SBC/sbc_oce.F90 (modified) (3 diffs)
• OPA_SRC/SBC/sbcmod.F90 (modified) (8 diffs)
• OPA_SRC/SBC/sbcrnf.F90 (modified) (8 diffs)
• OPA_SRC/TRA/trabbc.F90 (modified) (1 diff)
• OPA_SRC/TRA/trabbl.F90 (modified) (10 diffs)
• OPA_SRC/TRA/tranxt.F90 (modified) (17 diffs)
• OPA_SRC/TRA/traqsr.F90 (modified) (9 diffs)
• OPA_SRC/TRA/trasbc.F90 (modified) (8 diffs)
• OPA_SRC/istate.F90 (modified) (2 diffs)
• OPA_SRC/oce.F90 (modified) (1 diff)
• OPA_SRC/opa.F90 (modified) (1 diff)
• OPA_SRC/step.F90 (modified) (1 diff)
• OPA_SRC/step_oce.F90 (modified) (1 diff)
• TOP_SRC/PISCES/p4zmort.F90 (modified) (1 diff)
• TOP_SRC/PISCES/p4zrem.F90 (modified) (1 diff)
• TOP_SRC/PISCES/p4zsink.F90 (modified) (1 diff)
• TOP_SRC/TRP/trcadv.F90 (modified) (1 diff)
• TOP_SRC/TRP/trcdmp.F90 (modified) (1 diff)
• TOP_SRC/TRP/trcldf.F90 (modified) (1 diff)
• TOP_SRC/TRP/trcnxt.F90 (modified) (3 diffs)
• TOP_SRC/TRP/trcrad.F90 (modified) (1 diff)
• TOP_SRC/TRP/trcsbc.F90 (modified) (4 diffs)
• TOP_SRC/TRP/trczdf.F90 (modified) (1 diff)
• TOP_SRC/TRP/trdmld_trc.F90 (modified) (7 diffs)
• TOP_SRC/TRP/trdmod_trc.F90 (modified) (1 diff)
• TOP_SRC/oce_trc.F90 (modified) (1 diff)

branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DOM/dom_oce.F90

@@ -145 +145 @@
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3vw_1             !: analytical vertical scale factors at  VW--
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3w_1  , e3uw_1    !:                                       W--UW  points (m)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3t_b              !: before         -      -      -    -   T      points (m)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3u_b  , e3v_b     !:   -            -      -      -    -   U--V   points (m)
 #else
    LOGICAL, PUBLIC, PARAMETER ::   lk_vvl = .FALSE.   !: fixed grid flag

branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DOM/domvvl.F90

@@ -51 +51 @@
       !!----------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zcoefu, zcoefv, zcoeff
+      REAL(wp) ::   zcoefu , zcoefv   , zcoeff                   ! temporary scalars
+      REAL(wp) ::   zv_t_ij, zv_t_ip1j, zv_t_ijp1, zv_t_ip1jp1   !     -        -
+      REAL(wp), DIMENSION(jpi,jpj) ::  zs_t, zs_u_1, zs_v_1      !     -     2D workspace
       !!----------------------------------------------------------------------
 
@@ -62 +64 @@
       IF( lk_zco )   CALL ctl_stop( 'dom_vvl : key_zco is incompatible with variable volume option key_vvl')
 
-      IF( ln_zco) THEN
-         DO jk = 1, jpk
-            gdept(:,:,jk) = gdept_0(jk)
-            gdepw(:,:,jk) = gdepw_0(jk)
-            gdep3w(:,:,jk) = gdepw_0(jk)
-            e3t (:,:,jk) = e3t_0(jk)
-            e3u (:,:,jk) = e3t_0(jk)
-            e3v (:,:,jk) = e3t_0(jk)
-            e3f (:,:,jk) = e3t_0(jk)
-            e3w (:,:,jk) = e3w_0(jk)
-            e3uw(:,:,jk) = e3w_0(jk)
-            e3vw(:,:,jk) = e3w_0(jk)
-         END DO
-      ELSE
-         fsdept(:,:,:) = gdept (:,:,:)
-         fsdepw(:,:,:) = gdepw (:,:,:)
-         fsde3w(:,:,:) = gdep3w(:,:,:)
-         fse3t (:,:,:) = e3t   (:,:,:)
-         fse3u (:,:,:) = e3u   (:,:,:)
-         fse3v (:,:,:) = e3v   (:,:,:)
-         fse3f (:,:,:) = e3f   (:,:,:)
-         fse3w (:,:,:) = e3w   (:,:,:)
-         fse3uw(:,:,:) = e3uw  (:,:,:)
-         fse3vw(:,:,:) = e3vw  (:,:,:)
-      ENDIF
+      fsdept(:,:,:) = gdept (:,:,:)
+      fsdepw(:,:,:) = gdepw (:,:,:)
+      fsde3w(:,:,:) = gdep3w(:,:,:)
+      fse3t (:,:,:) = e3t   (:,:,:)
+      fse3u (:,:,:) = e3u   (:,:,:)
+      fse3v (:,:,:) = e3v   (:,:,:)
+      fse3f (:,:,:) = e3f   (:,:,:)
+      fse3w (:,:,:) = e3w   (:,:,:)
+      fse3uw(:,:,:) = e3uw  (:,:,:)
+      fse3vw(:,:,:) = e3vw  (:,:,:)
 
       !                                 !==  mu computation  ==!
@@ -130 +117 @@
          hv_0(:,:) = hv_0(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk)
       END DO
+      
+      ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
+      ! for ssh and scale factors
+      zs_t  (:,:) =       e1t(:,:) * e2t(:,:)
+      zs_u_1(:,:) = 0.5 / e1u(:,:) * e2u(:,:)
+      zs_v_1(:,:) = 0.5 / e1v(:,:) * e2v(:,:)
 
       DO jj = 1, jpjm1                          ! initialise before and now Sea Surface Height at u-, v-, f-points
          DO ji = 1, jpim1   ! NO vector opt.
-            zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
-            zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
-            zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
-            sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
-               &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )   
-            sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
-               &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )   
-            sshf_b(ji,jj) = zcoeff * ( sshb(ji  ,jj) + sshb(ji  ,jj+1)                 &
-               &                     + sshb(ji+1,jj) + sshb(ji+1,jj+1) )               
-            sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
-               &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )   
-            sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     & 
-               &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )     
-            sshf_n(ji,jj) = zcoeff * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)                 &
-               &                     + sshn(ji+1,jj) + sshn(ji+1,jj+1) )               
+            zcoefu = umask(ji,jj,1) * zs_u_1(ji,jj)
+            zcoefv = vmask(ji,jj,1) * zs_v_1(ji,jj)
+            zcoeff = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) / ( e1f(ji,jj) * e2f(ji,jj) )
+            ! before fields
+            zv_t_ij       = zs_t(ji  ,jj  ) * sshb(ji  ,jj  )
+            zv_t_ip1j     = zs_t(ji+1,jj  ) * sshb(ji+1,jj  )
+            zv_t_ijp1     = zs_t(ji  ,jj+1) * sshb(ji  ,jj+1)
+            sshu_b(ji,jj) = zcoefu * ( zv_t_ij + zv_t_ip1j )
+            sshv_b(ji,jj) = zcoefv * ( zv_t_ij + zv_t_ijp1 )
+            ! now fields
+            zv_t_ij       = zs_t(ji  ,jj  ) * sshn(ji  ,jj  )
+            zv_t_ip1j     = zs_t(ji+1,jj  ) * sshn(ji+1,jj  )
+            zv_t_ijp1     = zs_t(ji  ,jj+1) * sshn(ji  ,jj+1)
+            zv_t_ip1jp1   = zs_t(ji  ,jj+1) * sshn(ji  ,jj+1)
+            sshu_n(ji,jj) = zcoefu * ( zv_t_ij + zv_t_ip1j )
+            sshv_n(ji,jj) = zcoefv * ( zv_t_ij + zv_t_ijp1 )
+            sshf_n(ji,jj) = zcoeff * ( zv_t_ij + zv_t_ip1j + zv_t_ijp1 + zv_t_ip1jp1 )
          END DO
       END DO
-      CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )      ! lateral boundary conditions
-      CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
-      CALL lbc_lnk( sshf_b, 'F', 1. )   ;   CALL lbc_lnk( sshf_n, 'F', 1. )
+      CALL lbc_lnk( sshu_n, 'U', 1. )   ;   CALL lbc_lnk( sshu_b, 'U', 1. )      ! lateral boundary conditions
+      CALL lbc_lnk( sshv_n, 'V', 1. )   ;   CALL lbc_lnk( sshv_b, 'V', 1. )
+      CALL lbc_lnk( sshf_n, 'F', 1. )
+
+                                                ! initialise before scale factors at (u/v)-points
+      ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
+      DO jk = 1, jpkm1
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
+               zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
+               zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
+               fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) )
+               fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) )
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )               ! lateral boundary conditions
+      CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. )
+      ! Add initial scale factor to scale factor anomaly
+      fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
+      fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:)
       !
-         DO jk = 1, jpkm1
-            fsdept(:,:,jk) = fsdept_n(:,:,jk)          ! now local depths stored in fsdep. arrays
-            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
-            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
-            !
-            fse3t (:,:,jk) = fse3t_n (:,:,jk)          ! vertical scale factors stored in fse3. arrays
-            fse3u (:,:,jk) = fse3u_n (:,:,jk)
-            fse3v (:,:,jk) = fse3v_n (:,:,jk)
-            fse3f (:,:,jk) = fse3f_n (:,:,jk)
-            fse3w (:,:,jk) = fse3w_n (:,:,jk)
-            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
-            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
-         END DO
-
-
-
    END SUBROUTINE dom_vvl
 

branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DOM/domzgr_substitute.h90

@@ -46 +46 @@
 #   define  fse3vw(i,j,k)  e3vw_1(i,j,k)
 
-#   define  fsdept_b(i,j,k)  (fsdept_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k)))
-#   define  fsdepw_b(i,j,k)  (fsdepw_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k)))
-#   define  fsde3w_b(i,j,k)  (fsde3w_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k))-sshb(i,j))
-#   define  fse3t_b(i,j,k)   (fse3t_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k)))
-#   define  fse3u_b(i,j,k)   (fse3u_0(i,j,k)*(1.+sshu_b(i,j)*muu(i,j,k)))
-#   define  fse3v_b(i,j,k)   (fse3v_0(i,j,k)*(1.+sshv_b(i,j)*muv(i,j,k)))
-#   define  fse3f_b(i,j,k)   (fse3f_0(i,j,k)*(1.+sshf_b(i,j)*muf(i,j,k)))
-#   define  fse3w_b(i,j,k)   (fse3w_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k)))
+#   define  fse3t_b(i,j,k)   e3t_b(i,j,k)
+#   define  fse3u_b(i,j,k)   e3u_b(i,j,k)
+#   define  fse3v_b(i,j,k)   e3v_b(i,j,k)
 #   define  fse3uw_b(i,j,k)  (fse3uw_0(i,j,k)*(1.+sshu_b(i,j)*muu(i,j,k)))
 #   define  fse3vw_b(i,j,k)  (fse3vw_0(i,j,k)*(1.+sshv_b(i,j)*muv(i,j,k)))
@@ -70 +65 @@
 #   define  fse3t_m(i,j,k)   (fse3t_0(i,j,k)*(1.+ssh_m(i,j)*mut(i,j,k)))
 
-#   define  fsdept_a(i,j,k)  (fsdept_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k)))
-#   define  fsdepw_a(i,j,k)  (fsdepw_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k)))
-#   define  fsde3w_a(i,j,k)  (fsde3w_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k))-ssha(i,j))
 #   define  fse3t_a(i,j,k)   (fse3t_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k)))
 #   define  fse3u_a(i,j,k)   (fse3u_0(i,j,k)*(1.+sshu_a(i,j)*muu(i,j,k)))
 #   define  fse3v_a(i,j,k)   (fse3v_0(i,j,k)*(1.+sshv_a(i,j)*muv(i,j,k)))
-#   define  fse3f_a(i,j,k)   (fse3f_0(i,j,k)*(1.+sshf_a(i,j)*muf(i,j,k)))
-#   define  fse3w_a(i,j,k)   (fse3w_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k)))
-#   define  fse3uw_a(i,j,k)  (fse3uw_0(i,j,k)*(1.+sshu_a(i,j)*muu(i,j,k)))
-#   define  fse3vw_a(i,j,k)  (fse3vw_0(i,j,k)*(1.+sshv_a(i,j)*muv(i,j,k)))
 
 #else
@@ -94 +82 @@
 #   define  fse3vw(i,j,k)  fse3vw_0(i,j,k)
 
-#   define  fsdept_b(i,j,k)  fsdept_0(i,j,k)
-#   define  fsdepw_b(i,j,k)  fsdepw_0(i,j,k)
-#   define  fsde3w_b(i,j,k)  fsde3w_0(i,j,k)
 #   define  fse3t_b(i,j,k)   fse3t_0(i,j,k)
 #   define  fse3u_b(i,j,k)   fse3u_0(i,j,k)
 #   define  fse3v_b(i,j,k)   fse3v_0(i,j,k)
-#   define  fse3f_b(i,j,k)   fse3f_0(i,j,k)
-#   define  fse3w_b(i,j,k)   fse3w_0(i,j,k)
 #   define  fse3uw_b(i,j,k)  fse3uw_0(i,j,k)
 #   define  fse3vw_b(i,j,k)  fse3vw_0(i,j,k)
@@ -118 +101 @@
 #   define  fse3t_m(i,j,k)   fse3t_0(i,j,k)
 
-#   define  fsdept_a(i,j,k)  fsdept_0(i,j,k)
-#   define  fsdepw_a(i,j,k)  fsdepw_0(i,j,k)
-#   define  fsde3w_a(i,j,k)  fsde3w_0(i,j,k)
 #   define  fse3t_a(i,j,k)   fse3t_0(i,j,k)
 #   define  fse3u_a(i,j,k)   fse3u_0(i,j,k)
 #   define  fse3v_a(i,j,k)   fse3v_0(i,j,k)
-#   define  fse3f_a(i,j,k)   fse3f_0(i,j,k)
-#   define  fse3w_a(i,j,k)   fse3w_0(i,j,k)
-#   define  fse3uw_a(i,j,k)  fse3uw_0(i,j,k)
-#   define  fse3vw_a(i,j,k)  fse3vw_0(i,j,k)
 #endif
    !!----------------------------------------------------------------------

branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/divcur.F90

@@ -14 +14 @@
    USE in_out_manager  ! I/O manager
    USE obc_oce         ! ocean lateral open boundary condition
-   USE bdy_oce        ! Unstructured open boundaries variables
+   USE bdy_oce         ! Unstructured open boundaries variables
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
 
@@ -87 +87 @@
       INTEGER ::   ii, ij, jl     ! temporary integer
       INTEGER ::   ijt, iju       ! temporary integer
+      REAL(wp) ::  zraur, zdep
       REAL(wp), DIMENSION(   jpi  ,1:jpj+2) ::   zwu   ! workspace
       REAL(wp), DIMENSION(-1:jpi+2,  jpj  ) ::   zwv   ! workspace
@@ -301 +302 @@
       !! * Local declarations
       INTEGER  ::   ji, jj, jk          ! dummy loop indices
+      REAL(wp) ::  zraur, zdep
       !!----------------------------------------------------------------------
 

branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/dynnxt.F90

@@ -22 +22 @@
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
+   USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE phycst          ! physical constants
    USE dynspg_oce      ! type of surface pressure gradient
    USE dynadv          ! dynamics: vector invariant versus flux form
@@ -87 +89 @@
       !!               un,vn   now horizontal velocity of next time-step
       !!----------------------------------------------------------------------
+      USE oce, ONLY :   ze3u_f => ta   ! use ta as 3D workspace
+      USE oce, ONLY :   ze3v_f => sa   ! use sa as 3D workspace
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
       !!
@@ -95 +99 @@
       REAL(wp) ::   zue3a , zue3n , zue3b    ! temporary scalar
       REAL(wp) ::   zve3a , zve3n , zve3b    !    -         -
-      REAL(wp) ::   ze3u_b, ze3u_n, ze3u_a   !    -         -
-      REAL(wp) ::   ze3v_b, ze3v_n, ze3v_a   !    -         - 
       REAL(wp) ::   zuf   , zvf              !    -         - 
+      REAL(wp) ::   zec                      !    -         - 
+      REAL(wp) ::   zv_t_ij  , zv_t_ip1j     !     -        -
+      REAL(wp) ::   zv_t_ijp1                !     -        -
+      REAL(wp), DIMENSION(jpi,jpj) ::  zs_t, zs_u_1, zs_v_1      ! temporary 2D workspace
       !!----------------------------------------------------------------------
 
@@ -146 +152 @@
 # if defined key_obc
       !                                !* OBC open boundaries
-      IF( lk_obc )   CALL obc_dyn( kt )
+      CALL obc_dyn( kt )
       !
       IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN
@@ -212 +218 @@
          END DO
       ELSE                                             !* Leap-Frog : Asselin filter and swap
-         IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN          ! applied on velocity
+         !                                ! =============!
+         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
+            !                             ! =============!
             DO jk = 1, jpkm1                              
                DO jj = 1, jpj
                   DO ji = 1, jpi    
-                     zuf = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk)
-                     zvf = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk)
+                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
+                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
                      !
                      ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
@@ -226 +234 @@
                END DO
             END DO
-         ELSE                                                ! applied on thickness weighted velocity
+            !                             ! ================!
+         ELSE                             ! Variable volume !
+            !                             ! ================!
+            ! Before scale factor at t-points
+            ! -------------------------------
             DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     ze3u_a = fse3u_a(ji,jj,jk)
-                     ze3v_a = fse3v_a(ji,jj,jk)
-                     ze3u_n = fse3u_n(ji,jj,jk)
-                     ze3v_n = fse3v_n(ji,jj,jk)
-                     ze3u_b = fse3u_b(ji,jj,jk)
-                     ze3v_b = fse3v_b(ji,jj,jk)
-                     !
-                     zue3a = ua(ji,jj,jk) * ze3u_a
-                     zve3a = va(ji,jj,jk) * ze3v_a
-                     zue3n = un(ji,jj,jk) * ze3u_n
-                     zve3n = vn(ji,jj,jk) * ze3v_n
-                     zue3b = ub(ji,jj,jk) * ze3u_b
-                     zve3b = vb(ji,jj,jk) * ze3v_b
-                     !
-                     zuf  = ( atfp  * ( zue3b  + zue3a  ) + atfp1 * zue3n  )   &
-                        & / ( atfp  * ( ze3u_b + ze3u_a ) + atfp1 * ze3u_n ) * umask(ji,jj,jk)
-                     zvf  = ( atfp  * ( zve3b  + zve3a  ) + atfp1 * zve3n  )   &
-                        & / ( atfp  * ( ze3v_b + ze3v_a ) + atfp1 * ze3v_n ) * vmask(ji,jj,jk)
-                     !
-                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
-                     vb(ji,jj,jk) = zvf
-                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
-                     vn(ji,jj,jk) = va(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
+               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
+                  &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     &
+                  &                         - 2.e0 * fse3t_n(:,:,jk)            )
+            ENDDO
+            ! Add volume filter correction only at the first level of t-point scale factors
+            zec = atfp * rdt / rau0
+            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
+            ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
+            zs_t  (:,:) =       e1t(:,:) * e2t(:,:)
+            zs_u_1(:,:) = 0.5 / e1u(:,:) * e2u(:,:)
+            zs_v_1(:,:) = 0.5 / e1v(:,:) * e2v(:,:)
+            !
+            IF( ln_dynadv_vec ) THEN
+               ! Before scale factor at (u/v)-points
+               ! -----------------------------------
+               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
+               DO jk = 1, jpkm1
+                  DO jj = 1, jpjm1
+                     DO ji = 1, jpim1
+                        zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
+                        zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
+                        zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
+                        fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) )
+                        fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) )
+                     END DO
+                  END DO
+               END DO
+               ! lateral boundary conditions
+               CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )
+               CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. )
+               ! Add initial scale factor to scale factor anomaly
+               fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
+               fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:)
+               ! Leap-Frog - Asselin filter and swap: applied on velocity
+               ! -----------------------------------
+               DO jk = 1, jpkm1
+                  DO jj = 1, jpj
+                     DO ji = 1, jpi
+                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
+                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
+                        !
+                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
+                        vb(ji,jj,jk) = zvf
+                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
+                        vn(ji,jj,jk) = va(ji,jj,jk)
+                     END DO
+                  END DO
+               END DO
+               !
+            ELSE
+               ! Temporary filered scale factor at (u/v)-points (will become before scale factor)
+               !-----------------------------------------------
+               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
+               DO jk = 1, jpkm1
+                  DO jj = 1, jpjm1
+                     DO ji = 1, jpim1
+                        zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
+                        zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
+                        zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
+                        ze3u_f(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) )
+                        ze3v_f(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) )
+                     END DO
+                  END DO
+               END DO
+               ! lateral boundary conditions
+               CALL lbc_lnk( ze3u_f, 'U', 1. )
+               CALL lbc_lnk( ze3v_f, 'V', 1. )
+               ! Add initial scale factor to scale factor anomaly
+               ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:)
+               ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:)
+               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
+               ! -----------------------------------             ===========================
+               DO jk = 1, jpkm1
+                  DO jj = 1, jpj
+                     DO ji = 1, jpim1
+                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
+                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
+                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
+                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
+                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
+                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
+                        !
+                        zuf  = ( zue3n + atfp * ( zue3b - 2.e0 * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
+                        zvf  = ( zve3n + atfp * ( zve3b - 2.e0 * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
+                        !
+                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
+                        vb(ji,jj,jk) = zvf
+                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
+                        vn(ji,jj,jk) = va(ji,jj,jk)
+                     END DO
+                  END DO
+               END DO
+               fse3u_b(:,:,:) = ze3u_f(:,:,:)                   ! e3u_b <-- filtered scale factor
+               fse3v_b(:,:,:) = ze3v_f(:,:,:)
+               CALL lbc_lnk( ub, 'U', -1. )                     ! lateral boundary conditions
+               CALL lbc_lnk( vb, 'V', -1. )
+            ENDIF
+            !
          ENDIF
+         !
       ENDIF
 

branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/dynzdf_exp.F90

@@ -4 +4 @@
    !! Ocean dynamics:  vertical component(s) of the momentum mixing trend
    !!==============================================================================
-   !! History :      !  90-10  (B. Blanke)  Original code
-   !!                !  97-05  (G. Madec)  vertical component of isopycnal
-   !!           8.5  !  02-08  (G. Madec)  F90: Free form and module
+   !! History :  OPA  !  1990-10  (B. Blanke)  Original code
+   !!            8.0  !  1997-05  (G. Madec)  vertical component of isopycnal
+   !!   NEMO     1.0  !  1002-08  (G. Madec)  F90: Free form and module
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
    !!----------------------------------------------------------------------
 
@@ -13 +14 @@
    !!                  sion using an explicit time-stepping scheme.
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
@@ -24 +24 @@
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC dyn_zdf_exp    ! called by step.F90
+   PUBLIC   dyn_zdf_exp   ! called by step.F90
 
    !! * Substitutions
@@ -31 +30 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -47 +46 @@
       !!      technique. The vertical diffusion of momentum is given by:
       !!         diffu = dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ub) )
-      !!      Surface boundary conditions: wind stress input
+      !!      Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2)
       !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F90)
       !!      Add this trend to the general trend ua :
@@ -54 +53 @@
       !! ** Action : - Update (ua,va) with the vertical diffusive trend
       !!---------------------------------------------------------------------
-      !! * Arguments
-      INTEGER , INTENT( in ) ::   kt                           ! ocean time-step index
-      REAL(wp), INTENT( in ) ::   p2dt                         ! time-step 
-
-      !! * Local declarations
-      INTEGER ::   ji, jj, jk, jl                              ! dummy loop indices
-      REAL(wp) ::   zrau0r, zlavmr, zua, zva                   ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpk) ::   zwx, zwy, zwz, zww     ! temporary workspace arrays
+      INTEGER , INTENT(in) ::   kt     ! ocean time-step index
+      REAL(wp), INTENT(in) ::   p2dt   ! time-step 
+      !!
+      INTEGER ::   ji, jj, jk, jl                            ! dummy loop indices
+      REAL(wp) ::   zrau0r, zlavmr, zua, zva                 ! temporary scalars
+      REAL(wp), DIMENSION(jpi,jpk) ::   zwx, zwy, zwz, zww   ! 2D workspace
       !!----------------------------------------------------------------------
 
       IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn_zdf_exp : vertical momentum diffusion explicit operator'
+         IF(lwp) WRITE(numout,*) 'dyn_zdf_exp : vertical momentum diffusion - explicit operator'
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
       ENDIF
 
-      ! Local constant initialization
-      ! -----------------------------
-      zrau0r = 1. / rau0                                   ! inverse of the reference density
-      zlavmr = 1. / float( nn_zdfexp )                      ! inverse of the number of sub time step
+      zrau0r = 1. / rau0                              ! Local constant initialization
+      zlavmr = 1. / REAL( nn_zdfexp )
 
-      !                                                ! ===============
-      DO jj = 2, jpjm1                                 !  Vertical slab
-         !                                             ! ===============
-
-         ! Surface boundary condition
-         DO ji = 2, jpim1
-            zwy(ji,1) = utau(ji,jj) * zrau0r
-            zww(ji,1) = vtau(ji,jj) * zrau0r
+      !                                               ! ===============
+      DO jj = 2, jpjm1                                !  Vertical slab
+         !                                            ! ===============
+         DO ji = 2, jpim1         ! Surface boundary condition
+            zwy(ji,1) = ( utau_b(ji,jj) + utau(ji,jj) ) * zrau0r
+            zww(ji,1) = ( vtau_b(ji,jj) + vtau(ji,jj) ) * zrau0r
          END DO  
-
-         ! Initialization of x, z and contingently trends array
-         DO jk = 1, jpk
+         DO jk = 1, jpk         ! Initialization of x, z and contingently trends array
             DO ji = 2, jpim1
                zwx(ji,jk) = ub(ji,jj,jk)
@@ -92 +83 @@
             END DO  
          END DO  
-
-         ! Time splitting loop
-         DO jl = 1, nn_zdfexp
-
-            ! First vertical derivative
-            DO jk = 2, jpk
+         !
+         DO jl = 1, nn_zdfexp     ! Time splitting loop
+            !
+            DO jk = 2, jpk            ! First vertical derivative
                DO ji = 2, jpim1
                   zwy(ji,jk) = avmu(ji,jj,jk) * ( zwx(ji,jk-1) - zwx(ji,jk) ) / fse3uw(ji,jj,jk) 
@@ -103 +92 @@
                END DO  
             END DO  
-
-            ! Second vertical derivative and trend estimation at kt+l*rdt/nn_zdfexp
-            DO jk = 1, jpkm1
+            DO jk = 1, jpkm1          ! Second vertical derivative and trend estimation at kt+l*rdt/nn_zdfexp
                DO ji = 2, jpim1
-                  zua = zlavmr*( zwy(ji,jk) - zwy(ji,jk+1) ) / fse3u(ji,jj,jk)
-                  zva = zlavmr*( zww(ji,jk) - zww(ji,jk+1) ) / fse3v(ji,jj,jk)
+                  zua = zlavmr * ( zwy(ji,jk) - zwy(ji,jk+1) ) / fse3u(ji,jj,jk)
+                  zva = zlavmr * ( zww(ji,jk) - zww(ji,jk+1) ) / fse3v(ji,jj,jk)
                   ua(ji,jj,jk) = ua(ji,jj,jk) + zua
                   va(ji,jj,jk) = va(ji,jj,jk) + zva
-
-                  zwx(ji,jk) = zwx(ji,jk) + p2dt*zua*umask(ji,jj,jk)
-                  zwz(ji,jk) = zwz(ji,jk) + p2dt*zva*vmask(ji,jj,jk)
+                  !
+                  zwx(ji,jk) = zwx(ji,jk) + p2dt * zua * umask(ji,jj,jk)
+                  zwz(ji,jk) = zwz(ji,jk) + p2dt * zva * vmask(ji,jj,jk)
                END DO  
             END DO  
-
+            !
          END DO  
-
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
-
+         !                                            ! ===============
+      END DO                                          !   End of slab
+      !                                               ! ===============
    END SUBROUTINE dyn_zdf_exp
 

branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

@@ -4 +4 @@
    !! Ocean dynamics:  vertical component(s) of the momentum mixing trend
    !!==============================================================================
-   !! History :      !  90-10  (B. Blanke)  Original code
-   !!                !  97-05  (G. Madec)  vertical component of isopycnal
-   !!           8.5  !  02-08  (G. Madec)  F90: Free form and module
+   !! History :  OPA  !  1990-10  (B. Blanke)  Original code
+   !!            8.0  !  1997-05  (G. Madec)  vertical component of isopycnal
+   !!   NEMO     1.0  !  2002-08  (G. Madec)  F90: Free form and module
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
    !!----------------------------------------------------------------------
 
@@ -13 +14 @@
    !!                  sion using a implicit time-stepping.
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
-   !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
@@ -28 +24 @@
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC dyn_zdf_imp    ! called by step.F90
+   PUBLIC   dyn_zdf_imp   ! called by step.F90
 
    !! * Substitutions
@@ -35 +30 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
 CONTAINS
-
 
    SUBROUTINE dyn_zdf_imp( kt, p2dt )
@@ -54 +48 @@
       !!             dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
       !!      backward time stepping
-      !!      Surface boundary conditions: wind stress input
+      !!      Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2)
       !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F)
       !!      Add this trend to the general trend ua :
       !!         ua = ua + dz( avmu dz(u) )
       !!
-      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive
-      !!               mixing trend.
+      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend.
       !!---------------------------------------------------------------------
-      !! * Modules used
-      USE oce, ONLY :  zwd   => ta,   &                ! use ta as workspace
-                       zws   => sa                     ! use sa as workspace
-
-      !! * Arguments
-      INTEGER , INTENT( in ) ::   kt                   ! ocean time-step index
-      REAL(wp), INTENT( in ) ::  p2dt                  ! vertical profile of tracer time-step
-
-      !! * Local declarations
-      INTEGER ::   ji, jj, jk                          ! dummy loop indices
-      REAL(wp) ::   zrau0r, z2dtf, zcoef, zzws, zrhs   ! temporary scalars
-      REAL(wp) ::   zzwi                               ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi        ! temporary workspace arrays
+      USE oce, ONLY :  zwd   => ta      ! use ta as workspace
+      USE oce, ONLY :  zws   => sa      ! use sa as workspace
+      !!
+      INTEGER , INTENT( in ) ::   kt    ! ocean time-step index
+      REAL(wp), INTENT( in ) ::  p2dt   ! vertical profile of tracer time-step
+      !!
+      INTEGER  ::   ji, jj, jk             ! dummy loop indices
+      REAL(wp) ::   zrau0r, zcoef         ! temporary scalars
+      REAL(wp) ::   zzwi, zzws, zrhs       ! temporary scalars
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi        ! 3D workspace
       !!----------------------------------------------------------------------
 
@@ -95 +85 @@
       ! friction velocity in dyn_bfr using values from the previous timestep. There
       ! is no need to include these in the implicit calculation.
-      DO jk = 1, jpkm1
+      !
+      DO jk = 1, jpkm1        ! Matrix
          DO jj = 2, jpjm1 
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zcoef = - p2dt / fse3u(ji,jj,jk)
-               zzwi          = zcoef * avmu(ji,jj,jk  ) / fse3uw(ji,jj,jk  )
-               zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk)
-               zzws          = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
-               zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1)
+               zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
+               zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk)
+               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
+               zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1)
                zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws
             END DO
          END DO
       END DO
-
-      ! Surface boudary conditions
-      DO jj = 2, jpjm1 
+      DO jj = 2, jpjm1        ! Surface boudary conditions
          DO ji = fs_2, fs_jpim1   ! vector opt.
             zwi(ji,jj,1) = 0.
@@ -130 +119 @@
       !   The solution (the after velocity) is in ua
       !-----------------------------------------------------------------------
-
-      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
-      DO jk = 2, jpkm1
+      !
+      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -139 +127 @@
          END DO
       END DO
-
-      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
-      DO jj = 2, jpjm1   
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-!!! change les resultats (derniers digit, pas significativement + rapide 1* de moins)
-!!!         ua(ji,jj,1) = ub(ji,jj,1)  &
-!!!                      + p2dt * ( ua(ji,jj,1) + utau(ji,jj) / ( fse3u(ji,jj,1)*rau0 ) )
-            z2dtf = p2dt / ( fse3u(ji,jj,1)*rau0 )
-            ua(ji,jj,1) = ub(ji,jj,1)  &
-                         + p2dt *  ua(ji,jj,1) + z2dtf * utau(ji,jj)
+      !
+      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            ua(ji,jj,1) = ub(ji,jj,1)   &
+               &        + p2dt * (  ua(ji,jj,1) + 0.5 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( fse3u(ji,jj,1) * rau0 )  )
          END DO
       END DO
@@ -159 +142 @@
          END DO
       END DO
-
-      ! thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk
-      DO jj = 2, jpjm1   
+      !
+      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==
          DO ji = fs_2, fs_jpim1   ! vector opt.
             ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
@@ -173 +155 @@
          END DO
       END DO
-
-      ! Normalization to obtain the general momentum trend ua
-      DO jk = 1, jpkm1
+      !
+      DO jk = 1, jpkm1        !==  Normalization to obtain the general momentum trend ua  ==
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -192 +173 @@
       ! friction velocity in dyn_bfr using values from the previous timestep. There
       ! is no need to include these in the implicit calculation.
-      DO jk = 1, jpkm1
+      !
+      DO jk = 1, jpkm1        ! Matrix
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zcoef = -p2dt / fse3v(ji,jj,jk)
-               zzwi          = zcoef * avmv(ji,jj,jk  ) / fse3vw(ji,jj,jk  )
+               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  )
                zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk)
-               zzws       = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
+               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
                zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1)
                zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws
@@ -204 +186 @@
          END DO
       END DO
-
-      ! Surface boudary conditions
-      DO jj = 2, jpjm1   
+      DO jj = 2, jpjm1        ! Surface boudary conditions
          DO ji = fs_2, fs_jpim1   ! vector opt.
             zwi(ji,jj,1) = 0.e0
@@ -223 +203 @@
       !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
       !
-      !   m is decomposed in the product of an upper and lower triangular
-      !   matrix
+      !   m is decomposed in the product of an upper and lower triangular matrix
       !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
       !   The solution (after velocity) is in 2d array va
       !-----------------------------------------------------------------------
-
-      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
-      DO jk = 2, jpkm1
+      !
+      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -237 +215 @@
          END DO
       END DO
-
-      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-!!! change les resultats (derniers digit, pas significativement + rapide 1* de moins)
-!!!         va(ji,jj,1) = vb(ji,jj,1)  &
-!!!                      + p2dt * ( va(ji,jj,1) + vtau(ji,jj) / ( fse3v(ji,jj,1)*rau0 ) )
-            z2dtf = p2dt / ( fse3v(ji,jj,1)*rau0 )
-            va(ji,jj,1) = vb(ji,jj,1)  &
-                         + p2dt * va(ji,jj,1) + z2dtf * vtau(ji,jj)
+      !
+      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            va(ji,jj,1) = vb(ji,jj,1)   &
+               &        + p2dt * (  va(ji,jj,1) + 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( fse3v(ji,jj,1) * rau0 )  )
          END DO
       END DO
@@ -257 +230 @@
          END DO
       END DO
-
-      ! thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk
-      DO jj = 2, jpjm1   
+      !
+      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
          DO ji = fs_2, fs_jpim1   ! vector opt.
             va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
@@ -271 +243 @@
          END DO
       END DO
-
-      ! Normalization to obtain the general momentum trend va
-      DO jk = 1, jpkm1
+      !
+      DO jk = 1, jpkm1     !==  Normalization to obtain the general momentum trend va  ==
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -280 +251 @@
          END DO
       END DO
-
+      !
    END SUBROUTINE dyn_zdf_imp
 

branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/sshwzv.F90

@@ -5 +5 @@
    !!==============================================================================
    !! History :  3.1  !  2009-02  (G. Madec, M. Leclair)  Original code
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA 
    !!----------------------------------------------------------------------
 
@@ -27 +28 @@
    USE diaar5, ONLY :   lk_diaar5
    USE iom
-   USE sbcrnf, ONLY : rnf_dep, rnf_mod_dep  ! River runoff 
+   USE sbcrnf, ONLY : rnf_dep, rnf_mod_dep  ! River runoff
 
    IMPLICIT NONE
@@ -38 +39 @@
 #  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
-
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -54 +54 @@
       !!              and update the now vertical coordinate (lk_vvl=T).
       !!
-      !! ** Method  : -
-      !!
-      !!              - Using the incompressibility hypothesis, the vertical 
+      !! ** Method  : - Using the incompressibility hypothesis, the vertical 
       !!      velocity is computed by integrating the horizontal divergence  
       !!      from the bottom to the surface minus the scale factor evolution.
@@ -63 +61 @@
       !! ** action  :   ssha    : after sea surface height
       !!                wn      : now vertical velocity
-      !! if lk_vvl=T:   sshu_a, sshv_a, sshf_a  : after sea surface height
-      !!                          at u-, v-, f-point s
-      !!                hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points
+      !!                sshu_a, sshv_a, sshf_a  : after sea surface height (lk_vvl=T)
+      !!                hu, hv, hur, hvr        : ocean depth and its inverse at u-,v-points
+      !!
+      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
       !!----------------------------------------------------------------------
       USE oce, ONLY :   z3d => ta   ! use ta as 3D workspace
@@ -73 +72 @@
       INTEGER  ::   ji, jj, jk      ! dummy loop indices
       REAL(wp) ::   zcoefu, zcoefv, zcoeff      ! temporary scalars
-      REAL(wp) ::   z2dt, zraur     ! temporary scalars
+      REAL(wp) ::   z2dt, z1_2dt, z1_rau0       ! temporary scalars
       REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv       ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj) ::   z2d         ! 2D workspace
@@ -79 +78 @@
 
       IF( kt == nit000 ) THEN
+         !
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity '
@@ -95 +95 @@
                   sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
                      &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
-                  sshf_b(ji,jj) = zcoeff * ( sshb(ji  ,jj) + sshb(ji  ,jj+1)                 &
-                     &                     + sshb(ji+1,jj) + sshb(ji+1,jj+1) )
                   sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
                      &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
                   sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
                      &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
-                  sshf_n(ji,jj) = zcoeff * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)                 &
-                     &                     + sshn(ji+1,jj) + sshn(ji+1,jj+1) )
                END DO
             END DO
             CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
             CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
-            CALL lbc_lnk( sshf_b, 'F', 1. )   ;   CALL lbc_lnk( sshf_n, 'F', 1. )
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1      ! NO Vector Opt.
+                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
+                       &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
+                       &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
+                       &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
+               END DO
+            END DO
+            CALL lbc_lnk( sshf_n, 'F', 1. )
          ENDIF
          !
       ENDIF
 
-      !                                           !------------------------------!
-      IF( lk_vvl ) THEN                           !  Update Now Vertical coord.  !   (only in vvl case)
-      !                                           !------------------------------!
+      !                                           !------------------------------------------!
+      IF( lk_vvl ) THEN                           !  Regridding: Update Now Vertical coord.  !   (only in vvl case)
+         !                                        !------------------------------------------!
          DO jk = 1, jpkm1
-            fsdept(:,:,jk) = fsdept_n(:,:,jk)          ! now local depths stored in fsdep. arrays
+            fsdept(:,:,jk) = fsdept_n(:,:,jk)         ! now local depths stored in fsdep. arrays
             fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
             fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
             !
-            fse3t (:,:,jk) = fse3t_n (:,:,jk)          ! vertical scale factors stored in fse3. arrays
+            fse3t (:,:,jk) = fse3t_n (:,:,jk)         ! vertical scale factors stored in fse3. arrays
             fse3u (:,:,jk) = fse3u_n (:,:,jk)
             fse3v (:,:,jk) = fse3v_n (:,:,jk)
@@ -128 +132 @@
             fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
          END DO
-         !                                             ! now ocean depth (at u- and v-points)
-         hu(:,:) = hu_0(:,:) + sshu_n(:,:)
+         !
+         hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points)
          hv(:,:) = hv_0(:,:) + sshv_n(:,:)
-         !                                             ! now masked inverse of the ocean depth (at u- and v-points)
+         !                                            ! now masked inverse of the ocean depth (at u- and v-points)
          hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) )
          hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) )
          !
-         !
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               rnf_dep(ji,jj) = 0.
-               DO jk = 1, rnf_mod_dep(ji,jj)                          ! recalculates rnf_dep to be the depth 
-                  rnf_dep(ji,jj) = rnf_dep(ji,jj) + fse3t(ji,jj,jk)    ! in metres to the bottom of the relevant grid box 
-               ENDDO
-            ENDDO
-         ENDDO
-         ! 
-      ENDIF
-
-                         CALL div_cur( kt )            ! Horizontal divergence & Relative vorticity
-      IF( n_cla == 1 )   CALL div_cla( kt )            ! Cross Land Advection (Update Hor. divergence)
-
-      ! set time step size (Euler/Leapfrog)
-      z2dt = 2. * rdt 
+      ENDIF
+      !
+
+                         CALL div_cur( kt )           ! Horizontal divergence & Relative vorticity
+      IF( n_cla == 1 )   CALL div_cla( kt )           ! Cross Land Advection (Update Hor. divergence)
+
+      z2dt = 2. * rdt                                 ! set time step size (Euler/Leapfrog)
       IF( neuler == 0 .AND. kt == nit000 )   z2dt =rdt
-
-      zraur = 1. / rau0
 
       !                                           !------------------------------!
@@ -163 +155 @@
         zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
       END DO
-
       !                                                ! Sea surface elevation time stepping
-      ssha(:,:) = (  sshb(:,:) - z2dt * ( zraur * ( emp(:,:)-rnf(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1) 
+      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
+      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b -2*rnf ) = emp
+      z1_rau0 = 0.5 / rau0
+      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) - 2 * rnf(:,:) ) + zhdiv(:,:) )  ) &
+      &                      * tmask(:,:,1)
 
 #if defined key_obc
-      IF ( Agrif_Root() ) THEN 
+      IF( Agrif_Root() ) THEN 
          ssha(:,:) = ssha(:,:) * obctmsk(:,:)
-         CALL lbc_lnk( ssha, 'T', 1. )  ! absolutly compulsory !! (jmm)
+         CALL lbc_lnk( ssha, 'T', 1. )                ! absolutly compulsory !! (jmm)
       ENDIF
 #endif
-
       !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
       IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
@@ -184 +178 @@
                   &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
                   &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
-               sshf_a(ji,jj) = 0.25 * umask(ji,jj,1) * umask (ji,jj+1,1)                                 & 
-                  &                                  * ( ssha(ji  ,jj) + ssha(ji  ,jj+1)                 &
-                  &                                    + ssha(ji+1,jj) + ssha(ji+1,jj+1) )
             END DO
          END DO
-         CALL lbc_lnk( sshu_a, 'U', 1. )               ! Boundaries conditions
+         ! Boundaries conditions
+         CALL lbc_lnk( sshu_a, 'U', 1. )
          CALL lbc_lnk( sshv_a, 'V', 1. )
-         CALL lbc_lnk( sshf_a, 'F', 1. )
-      ENDIF
-
+      ENDIF
       !                                           !------------------------------!
       !                                           !     Now Vertical Velocity    !
       !                                           !------------------------------!
-      !                                                ! integrate from the bottom the hor. divergence
-      DO jk = jpkm1, 1, -1
-         wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)   &
-              &                    - (  fse3t_a(:,:,jk)                   &
-              &                       - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt
+      z1_2dt = 1.e0 / z2dt
+      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
+         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
+         wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
+            &                      - (  fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
+            &                         * tmask(:,:,jk) * z1_2dt
       END DO
-      !
+
+      !                                           !------------------------------!
+      !                                           !           outputs            !
+      !                                           !------------------------------!
       CALL iom_put( "woce", wn                    )   ! vertical velocity
       CALL iom_put( "ssh" , sshn                  )   ! sea surface height
       CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
-      IF( lk_diaar5 ) THEN
+      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
+         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
          z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
          DO jk = 1, jpk
             z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
          END DO
-         CALL iom_put( "w_masstr" , z3d                     )   !           vertical mass transport
-         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )   ! square of vertical mass transport
-      ENDIF
+         CALL iom_put( "w_masstr" , z3d                     )  
+         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
+      ENDIF
+      !
+      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
       !
    END SUBROUTINE ssh_wzv
@@ -227 +224 @@
       !!              ssha  already computed in ssh_wzv  
       !!
-      !! ** Method  : - apply Asselin time fiter to now ssh and swap :
-      !!             sshn = ssha + atfp * ( sshb -2 sshn + ssha )
-      !!             sshn = ssha
+      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
+      !!              from the filter, see Leclair and Madec 2010) and swap :
+      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
+      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
+      !!                sshn = ssha
       !!
       !! ** action  : - sshb, sshn   : before & now sea surface height
       !!                               ready for the next time step
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
-      !!
-      INTEGER  ::   ji, jj               ! dummy loop indices
+      !!
+      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !!
+      INTEGER  ::   ji, jj   ! dummy loop indices
+      REAL(wp) ::   zec      ! temporary scalar
       !!----------------------------------------------------------------------
 
@@ -245 +247 @@
       ENDIF
 
-      ! Time filter and swap of the ssh
-      ! -------------------------------
-
-      IF( lk_vvl ) THEN      ! Variable volume levels :   ssh at t-, u-, v, f-points
-         !                   ! ---------------------- !
+      !                       !--------------------------!
+      IF( lk_vvl ) THEN       !  Variable volume levels  !
+         !                    !--------------------------!
+         !
+         ! ssh at t-, u-, v, f-points
+         !===========================
          IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
             sshn  (:,:) = ssha  (:,:)                        ! now <-- after  (before already = now)
             sshu_n(:,:) = sshu_a(:,:)
             sshv_n(:,:) = sshv_a(:,:)
-            sshf_n(:,:) = sshf_a(:,:)
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1      ! NO Vector Opt.
+                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
+                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
+                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
+                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
+               END DO
+            END DO
+            ! Boundaries conditions
+            CALL lbc_lnk( sshf_n, 'F', 1. )
          ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
+            zec = atfp * rdt / rau0
             DO jj = 1, jpj
                DO ji = 1, jpi                                ! before <-- now filtered
-                  sshb  (ji,jj) = sshn(ji,jj)   + atfp * ( sshb  (ji,jj) - 2 * sshn  (ji,jj) + ssha  (ji,jj) )
-                  sshu_b(ji,jj) = sshu_n(ji,jj) + atfp * ( sshu_b(ji,jj) - 2 * sshu_n(ji,jj) + sshu_a(ji,jj) )
-                  sshv_b(ji,jj) = sshv_n(ji,jj) + atfp * ( sshv_b(ji,jj) - 2 * sshv_n(ji,jj) + sshv_a(ji,jj) )
-                  sshf_b(ji,jj) = sshf_n(ji,jj) + atfp * ( sshf_b(ji,jj) - 2 * sshf_n(ji,jj) + sshf_a(ji,jj) )
+                  sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   &
+                     &                          - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1)
                   sshn  (ji,jj) = ssha  (ji,jj)              ! now <-- after
                   sshu_n(ji,jj) = sshu_a(ji,jj)
                   sshv_n(ji,jj) = sshv_a(ji,jj)
-                  sshf_n(ji,jj) = sshf_a(ji,jj)
-               END DO
-            END DO
+               END DO
+            END DO
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1      ! NO Vector Opt.
+                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
+                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
+                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
+                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
+               END DO
+            END DO
+            ! Boundaries conditions
+            CALL lbc_lnk( sshf_n, 'F', 1. )
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1      ! NO Vector Opt.
+                  sshu_b(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
+                     &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
+                     &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
+                  sshv_b(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
+                     &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
+                     &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
+               END DO
+            END DO
+            ! Boundaries conditions
+            CALL lbc_lnk( sshu_b, 'U', 1. )
+            CALL lbc_lnk( sshv_b, 'V', 1. )
          ENDIF
-         !
-      ELSE                   ! fixed levels :   ssh at t-point only
-         !                   ! ------------ !
+         !                    !--------------------------!
+      ELSE                    !        fixed levels      !
+         !                    !--------------------------!
+         !
+         ! ssh at t-point only
+         !====================
          IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
             sshn(:,:) = ssha(:,:)                            ! now <-- after  (before already = now)
@@ -278 +314 @@
             DO jj = 1, jpj
                DO ji = 1, jpi                                ! before <-- now filtered
-                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )    
+                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )
                   sshn(ji,jj) = ssha(ji,jj)                  ! now <-- after
                END DO
@@ -286 +322 @@
       ENDIF
       !
-      IF(ln_ctl)   CALL prt_ctl(tab2d_1=sshb    , clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
+      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
       !
    END SUBROUTINE ssh_nxt

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

@@ -7 +7 @@
    !!   NEMO     1.0  !  2002-08  (G. Madec)  F90: Free form
    !!            2.0  !  2006-07  (S. Masson)  use IOM for restart
-   !!            3.3  !  2010-10  (C. Ethe, G. Madec) TRC-TRA merge (T-S in 4D)
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA
+   !!            - -  !  2010-10  (C. Ethe, G. Madec) TRC-TRA merge (T-S in 4D)
    !!----------------------------------------------------------------------
 
@@ -26 +27 @@
    USE zdfmxl          ! mixed layer depth
    USE trdmld_oce      ! ocean active mixed layer tracers trends variables
+   USE domvvl          ! variable volume
    USE traswp          ! swap from 4D T-S to 3D T & S and vice versa
 
@@ -39 +41 @@
 
    !! * Substitutions
+#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -123 +126 @@
       CALL iom_rstput( kt, nitrst, numrow, 'hdivb'  , hdivb   )
       CALL iom_rstput( kt, nitrst, numrow, 'sshb'   , sshb    )
+      IF( lk_vvl ) &
+      CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) )
       !
       CALL iom_rstput( kt, nitrst, numrow, 'un'     , un      )     ! now fields
@@ -154 +159 @@
       !!----------------------------------------------------------------------
       REAL(wp) ::   zrdt, zrdttra1
-      INTEGER  ::   jlibalt = jprstlib
+      INTEGER  ::   jk, jlibalt = jprstlib
       LOGICAL  ::   llok
       !!----------------------------------------------------------------------
@@ -192 +197 @@
       CALL iom_get( numror, jpdom_autoglo, 'hdivb', hdivb )
       CALL iom_get( numror, jpdom_autoglo, 'sshb' , sshb  )
+      IF( lk_vvl )  &
+      CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
       !
       CALL iom_get( numror, jpdom_autoglo, 'un'   , un    )        ! now    fields
@@ -219 +226 @@
          hdivb(:,:,:) = hdivn(:,:,:)
          sshb (:,:)   = sshn (:,:)
+         IF( lk_vvl ) THEN
+            DO jk = 1, jpk
+               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
+            ENDDO
+         ENDIF
       ENDIF
       !

branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/SBC/sbc_oce.F90

@@ -6 +6 @@
    !! History :  3.0   !  2006-06  (G. Madec)  Original code
    !!             -    !  2008-08  (G. Madec)  namsbc moved from sbcmod
+   !!            3.3   !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
    !!----------------------------------------------------------------------
    USE par_oce          ! ocean parameters
@@ -37 +38 @@
    !!              Ocean Surface Boundary Condition fields
    !!----------------------------------------------------------------------
-   LOGICAL , PUBLIC ::   lhftau = .FALSE.              !: HF tau contribution: mean of stress module - module of the mean stress
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   utau      !: sea surface i-stress (ocean referential)     [N/m2]
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   vtau      !: sea surface j-stress (ocean referential)     [N/m2]
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   taum      !: module of sea surface stress (at T-point)    [N/m2] 
-   !! wndm is used only in PISCES to compute gases exchanges at the surface of the free ocean or in the leads in sea-ice parts
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   wndm      !: wind speed module at T-point (=|U10m-Uoce|)  [m/s] 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qsr       !: sea heat flux:     solar                     [W/m2]
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qns       !: sea heat flux: non solar                     [W/m2]
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qsr_tot   !: total     solar heat flux (over sea and ice) [W/m2]
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qns_tot   !: total non solar heat flux (over sea and ice) [W/m2]
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp       !: freshwater budget: volume flux               [Kg/m2/s]
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emps      !: freshwater budget: concentration/dillution   [Kg/m2/s]
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rnf       !: river runoff   [Kg/m2/s]  
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp_tot   !: total evaporation - (liquid + solid) precpitation over oce and ice
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   tprecip   !: total precipitation           [Kg/m2/s]
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sprecip   !: solid precipitation           [Kg/m2/s]
-!!$   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rrunoff       !: runoff
-!!$   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   calving       !: calving
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fr_i      !: ice fraction  (between 0 to 1)               -
+   LOGICAL , PUBLIC ::   lhftau = .FALSE.              !: HF tau used in TKE: mean(stress module) - module(mean stress)
+   !!                                   !!   now    ! before   !!
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   utau   , utau_b   !: sea surface i-stress (ocean referential)     [N/m2]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   vtau   , vtau_b   !: sea surface j-stress (ocean referential)     [N/m2]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   taum              !: module of sea surface stress (at T-point)    [N/m2] 
+   !! wndm is used only in PISCES to compute surface gases exchanges in ice-free ocean or leads
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   wndm              !: wind speed module at T-point (=|U10m-Uoce|)  [m/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qsr               !: sea heat flux:     solar                     [W/m2]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qns    , qns_b    !: sea heat flux: non solar                     [W/m2]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qsr_tot           !: total     solar heat flux (over sea and ice) [W/m2]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qns_tot           !: total non solar heat flux (over sea and ice) [W/m2]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp    , emp_b    !: freshwater budget: volume flux               [Kg/m2/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emps   , emps_b   !: freshwater budget: concentration/dillution   [Kg/m2/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp_tot           !: total E-P over ocean and ice                 [Kg/m2/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rnf               !: river runoff   [Kg/m2/s]  
+   ! - ML - begin
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sbc_hc_n          !: sbc heat content trend now                   [K.m/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sbc_hc_b          !:  "   "      "      "   before                   "
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sbc_sc_n          !: sbc salt content trend now                   [psu.m/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sbc_sc_b          !:  "   "      "      "   before                   "
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   qsr_hc_n      !: heat content trend due to qsr flux now       [K.m/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   qsr_hc_b      !:  "      "      "    "  "   "   "   before       "
+   ! - ML - end
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   tprecip           !: total precipitation                          [Kg/m2/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sprecip           !: solid precipitation                          [Kg/m2/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fr_i              !: ice fraction = 1 - lead fraction      (between 0 to 1)
 #if defined key_cpl_carbon_cycle
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   atm_co2   !: atmospheric pCO2                             [ppm]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   atm_co2           !: atmospheric pCO2                             [ppm]
 #endif
+!!$   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rrunoff        !: runoff
+!!$   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   calving        !: calving
 
    !!----------------------------------------------------------------------
@@ -71 +81 @@
 
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2006) 
-   !! $ Id: $
+   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 
+   !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!======================================================================
-
 END MODULE sbc_oce

branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/SBC/sbcmod.F90

@@ -4 +4 @@
    !! Surface module :  provide to the ocean its surface boundary condition
    !!======================================================================
-   !! History :  3.0   !  07-2006  (G. Madec)  Original code
-   !!             -    !  08-2008  (S. Masson, E. .... ) coupled interface
+   !! History :  3.0  !  2006-07  (G. Madec)  Original code
+   !!            3.1  !  2008-08  (S. Masson, E. Maisonnave, G. Madec) coupled interface
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
    !!----------------------------------------------------------------------
 
@@ -49 +50 @@
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 
+   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -86 +87 @@
 !!gm  IF( lk_sbc_cpl       ) THEN   ;   ln_cpl      = .TRUE.   ;   ELSE   ;   ln_cpl      = .FALSE.   ;   ENDIF
 
-      IF ( Agrif_Root() ) THEN
+      IF( Agrif_Root() ) THEN
         IF( lk_lim2 )            nn_ice      = 2
         IF( lk_lim3 )            nn_ice      = 3
@@ -179 +180 @@
       !!                CAUTION : never mask the surface stress field (tke sbc)
       !!
-      !! ** Action  : - set the ocean surface boundary condition, i.e.  
-      !!                utau, vtau, qns, qsr, emp, emps, qrp, erp
+      !! ** Action  : - set the ocean surface boundary condition at before and now 
+      !!                time step, i.e.  
+      !!                utau_b, vtau_b, qns_b, qsr_b, emp_n, emps_b, qrp_b, erp_b
+      !!                utau  , vtau  , qns  , qsr  , emp  , emps  , qrp  , erp
       !!              - updte the ice fraction : fr_i
       !!----------------------------------------------------------------------
@@ -186 +189 @@
       !!---------------------------------------------------------------------
 
-      CALL iom_setkt( kt + nn_fsbc - 1 )         !  in sbc, iom_put is called every nn_fsbc time step
-      !
-      ! ocean to sbc mean sea surface variables (ss._m)
-      ! ---------------------------------------
-      CALL sbc_ssm( kt )                         ! sea surface mean currents (at U- and V-points), 
-      !                                          ! temperature and salinity (at T-point) over nf_sbc time-step
-      !                                          ! (i.e. sst_m, sss_m, ssu_m, ssv_m)
-
-      ! sbc formulation
-      ! ---------------
-         
-      SELECT CASE( nsbc )                        ! Compute ocean surface boundary condition
-      !                                          ! (i.e. utau,vtau, qns, qsr, emp, emps)
+      !                                            ! ---------------------------------------- !
+      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
+         !                                         ! ---------------------------------------- !
+         utau_b(:,:) = utau(:,:)                         ! Swap the ocean forcing fields
+         vtau_b(:,:) = vtau(:,:)                         ! (except at nit000 where before fields
+         qns_b (:,:) = qns (:,:)                         !  are set at the end of the routine)
+         ! The 3D heat content due to qsr forcing is treated in traqsr
+         ! qsr_b (:,:) = qsr (:,:)
+         emp_b (:,:) = emp (:,:)
+         emps_b(:,:) = emps(:,:)
+      ENDIF
+      !                                            ! ---------------------------------------- !
+      !                                            !        forcing field computation         !
+      !                                            ! ---------------------------------------- !
+
+      CALL iom_setkt( kt + nn_fsbc - 1 )                 ! in sbc, iom_put is called every nn_fsbc time step
+      !
+      CALL sbc_ssm( kt )                                 ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m)
+      !                                                  ! averaged over nf_sbc time-step
+
+                                                   !==  sbc formulation  ==!
+                                                            
+      SELECT CASE( nsbc )                                ! Compute ocean surface boundary condition
+      !                                                  ! (i.e. utau,vtau, qns, qsr, emp, emps)
       CASE(  0 )   ;   CALL sbc_gyre    ( kt )                    ! analytical formulation : GYRE configuration
       CASE(  1 )   ;   CALL sbc_ana     ( kt )                    ! analytical formulation : uniform sbc
@@ -214 +228 @@
       END SELECT
 
-      ! Misc. Options
-      ! -------------
+      !                                            !==  Misc. Options  ==!
 
 !!gm  IF( ln_dm2dc       )   CALL sbc_dcy( kt )                 ! Daily mean qsr distributed over the Diurnal Cycle
@@ -236 +249 @@
       !                                                         ! (update freshwater fluxes)
       !
+      IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
+         !                                             ! ---------------------------------------- !
+         IF( ln_rstart .AND.    &                               !* Restart: read in restart file
+            & iom_varid( numror, 'utau_b', ldstop = .FALSE. ) > 0 ) THEN 
+            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields red in the restart file'
+            CALL iom_get( numror, jpdom_autoglo, 'utau_b', utau_b )   ! before i-stress  (U-point)
+            CALL iom_get( numror, jpdom_autoglo, 'vtau_b', vtau_b )   ! before j-stress  (V-point)
+            CALL iom_get( numror, jpdom_autoglo, 'qns_b' , qns_b  )   ! before non solar heat flux (T-point)
+            ! The 3D heat content due to qsr forcing is treated in traqsr
+            ! CALL iom_get( numror, jpdom_autoglo, 'qsr_b' , qsr_b  )   ! before     solar heat flux (T-point)
+            CALL iom_get( numror, jpdom_autoglo, 'emp_b' , emp_b  )   ! before     freshwater flux (T-point)
+            CALL iom_get( numror, jpdom_autoglo, 'emps_b', emps_b )   ! before C/D freshwater flux (T-point)
+         ELSE                                                   !* no restart: set from nit000 values
+            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields set to nit000'
+            utau_b(:,:) = utau(:,:) 
+            vtau_b(:,:) = vtau(:,:)
+            qns_b (:,:) = qns (:,:)
+            ! qsr_b (:,:) = qsr (:,:)
+            emp_b (:,:) = emp (:,:)
+            emps_b(:,:) = emps(:,:)
+         ENDIF
+      ENDIF
+      !                                                ! ---------------------------------------- !
+      IF( lrst_oce ) THEN                              !      Write in the ocean restart file     !
+         !                                             ! ---------------------------------------- !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'sbc : ocean surface forcing fields written in ocean restart file ',   &
+            &                    'at it= ', kt,' date= ', ndastp
+         IF(lwp) WRITE(numout,*) '~~~~'
+         CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau )
+         CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtau )
+         CALL iom_rstput( kt, nitrst, numrow, 'qns_b'  , qns  )
+         ! The 3D heat content due to qsr forcing is treated in traqsr
+         ! CALL iom_rstput( kt, nitrst, numrow, 'qsr_b'  , qsr  )
+         CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp  )
+         CALL iom_rstput( kt, nitrst, numrow, 'emps_b' , emps )
+      ENDIF
+
+      !                                                ! ---------------------------------------- !
+      !                                                !        Outputs and control print         !
+      !                                                ! ---------------------------------------- !
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
-         CALL iom_put( "emp-rnf"  , (emp-rnf)  )                ! upward water flux 
-         CALL iom_put( "emps-rnf" , (emps-rnf) )                ! c/d water flux 
-         CALL iom_put( "qns+qsr"  , qns + qsr  )                ! total heat flux   (caution if ln_dm2dc=true, to be 
-         CALL iom_put( "qns"      , qns        )                ! solar heat flux    moved after the call to iom_setkt)
-         CALL iom_put( "qsr"      ,       qsr  )                ! solar heat flux    moved after the call to iom_setkt)
-         IF(  nn_ice > 0 )   CALL iom_put( "ice_cover", fr_i )  ! ice fraction 
+         CALL iom_put( "emp-rnf" , emp  - rnf )                   ! upward water flux
+         CALL iom_put( "emps-rnf", emps - rnf )                   ! c/d water flux
+         CALL iom_put( "qns+qsr" , qns  + qsr )                   ! total heat flux   (caution if ln_dm2dc=true, to be 
+         CALL iom_put( "qns"     , qns        )                   ! solar heat flux    moved after the call to iom_setkt)
+         CALL iom_put( "qsr"     ,       qsr  )                   ! solar heat flux    moved after the call to iom_setkt)
+         IF( nn_ice > 0 )   CALL iom_put( "ice_cover", fr_i )   ! ice fraction 
       ENDIF
       !
@@ -253 +307 @@
       !
       IF(ln_ctl) THEN         ! print mean trends (used for debugging)
-         CALL prt_ctl(tab2d_1=fr_i   , clinfo1=' fr_i - : ', mask1=tmask, ovlap=1 )
-         CALL prt_ctl(tab2d_1=(emp-rnf) , clinfo1=' emp-rnf  - : ', mask1=tmask, ovlap=1 ) 
-         CALL prt_ctl(tab2d_1=(emps-rnf), clinfo1=' emps-rnf - : ', mask1=tmask, ovlap=1 ) 
-         CALL prt_ctl(tab2d_1=qns    , clinfo1=' qns  - : ', mask1=tmask, ovlap=1 )
-         CALL prt_ctl(tab2d_1=qsr    , clinfo1=' qsr  - : ', mask1=tmask, ovlap=1 )
-         CALL prt_ctl(tab3d_1=tmask  , clinfo1=' tmask  : ', mask1=tmask, ovlap=1, kdim=jpk )
-         CALL prt_ctl(tab3d_1=tn     , clinfo1=' sst  - : ', mask1=tmask, ovlap=1, kdim=1   )
-         CALL prt_ctl(tab3d_1=sn     , clinfo1=' sss  - : ', mask1=tmask, ovlap=1, kdim=1   )
-         CALL prt_ctl(tab2d_1=utau   , clinfo1=' utau - : ', mask1=umask,                      &
-            &         tab2d_2=vtau   , clinfo2=' vtau - : ', mask2=vmask, ovlap=1 )
+         CALL prt_ctl(tab2d_1=fr_i      , clinfo1=' fr_i     - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab2d_1=(emp-rnf) , clinfo1=' emp-rnf  - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab2d_1=(emps-rnf), clinfo1=' emps-rnf - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab2d_1=qns       , clinfo1=' qns      - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab2d_1=qsr       , clinfo1=' qsr      - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab3d_1=tmask     , clinfo1=' tmask    - : ', mask1=tmask, ovlap=1, kdim=jpk )
+         CALL prt_ctl(tab3d_1=tn        , clinfo1=' sst      - : ', mask1=tmask, ovlap=1, kdim=1   )
+         CALL prt_ctl(tab3d_1=sn        , clinfo1=' sss      - : ', mask1=tmask, ovlap=1, kdim=1   )
+         CALL prt_ctl(tab2d_1=utau      , clinfo1=' utau     - : ', mask1=umask,                      &
+            &         tab2d_2=vtau      , clinfo2=' vtau     - : ', mask2=vmask, ovlap=1 )
       ENDIF
       !

branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/SBC/sbcrnf.F90

@@ -33 +33 @@
    TYPE(FLD_N)       , PUBLIC ::   sn_cnf                 !: information about the runoff mouth file to be read
    TYPE(FLD_N)                ::   sn_sal_rnf             !: information about the salinities of runoff file to be read  
-   TYPE(FLD_N)                ::   sn_tem_rnf             !: information about the temperatures of runoff file to be read  
+   TYPE(FLD_N)                ::   sn_tmp_rnf             !: information about the temperatures of runoff file to be read  
    TYPE(FLD_N)                ::   sn_dep_rnf             !: information about the depth which river inflow affects
    LOGICAL           , PUBLIC ::   ln_rnf_mouth = .false. !: specific treatment in mouths vicinity
@@ -53 +53 @@
    INTEGER, PUBLIC, DIMENSION(jpi,jpj) ::  rnf_mod_dep     !: depth of runoff in model levels  
    REAL,    PUBLIC, DIMENSION(jpi,jpj) ::  rnf_sal         !: salinity of river runoff  
-   REAL,    PUBLIC, DIMENSION(jpi,jpj) ::  rnf_tem         !: temperature of river runoff  
+   REAL,    PUBLIC, DIMENSION(jpi,jpj) ::  rnf_tmp         !: temperature of river runoff  
   
    INTEGER  ::  ji, jj ,jk    ! dummy loop indices  
@@ -86 +86 @@
       !!----------------------------------------------------------------------
       !                                   
-      IF( kt == nit000 ) THEN  
-         CALL sbc_rnf_init                      ! Read namelist and allocate structures
-      ENDIF
+      IF( kt == nit000 )  CALL sbc_rnf_init     ! Read namelist and allocate structures 
 
       !                                                   !-------------------!
@@ -117 +115 @@
             IF ( ln_rnf_att ) THEN  
                rnf_sal(:,:) = ( sf_sal_rnf(1)%fnow(:,:,1) )  
-               rnf_tem(:,:) = ( sf_tem_rnf(1)%fnow(:,:,1) )  
+               rnf_tmp(:,:) = ( sf_tem_rnf(1)%fnow(:,:,1) )  
             ELSE  
                rnf_sal(:,:) = 0  
-               rnf_tem(:,:) = -999  
+               rnf_tmp(:,:) = -999  
             ENDIF  
             CALL iom_put( "runoffs", rnf )         ! runoffs
@@ -143 +141 @@
       CHARACTER(len=32) ::   rn_dep_file   ! runoff file name  
       !! 
-      NAMELIST/namsbc_rnf/ cn_dir, ln_rnf_emp, sn_rnf, sn_cnf, sn_sal_rnf, sn_tem_rnf, sn_dep_rnf,   &  
+      NAMELIST/namsbc_rnf/ cn_dir, ln_rnf_emp, sn_rnf, sn_cnf, sn_sal_rnf, sn_tmp_rnf, sn_dep_rnf,   &  
          &                 ln_rnf_mouth, ln_rnf_att, rn_hrnf, rn_avt_rnf, rn_rfact  
       !!----------------------------------------------------------------------
@@ -157 +155 @@
 
       sn_sal_rnf = FLD_N( 'runoffs',  24.  , 'rosaline' ,  .TRUE.    , .true. ,   'yearly'  , ''    , ''  )  
-      sn_tem_rnf = FLD_N( 'runoffs',  24.  , 'rotemper' ,  .TRUE.    , .true. ,   'yearly'  , ''    , ''  )  
+      sn_tmp_rnf = FLD_N( 'runoffs',  24.  , 'rotemper' ,  .TRUE.    , .true. ,   'yearly'  , ''    , ''  )  
       sn_dep_rnf = FLD_N( 'runoffs',   0.  , 'rodepth'  ,  .FALSE.   , .true. ,   'yearly'  , ''    , ''  )  
       !
@@ -218 +216 @@
          IF ( ln_rnf_att ) THEN  
             CALL fld_fill (sf_sal_rnf, (/ sn_sal_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff salinity data', 'namsbc_rnf' )  
-            CALL fld_fill (sf_tem_rnf, (/ sn_tem_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff temperature data', 'namsbc_rnf' )  
+            CALL fld_fill (sf_tem_rnf, (/ sn_tmp_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff temperature data', 'namsbc_rnf' )  
   
             rn_dep_file = TRIM( cn_dir )//TRIM( sn_dep_rnf%clname )  
@@ -248 +246 @@
       ENDIF
       !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            rnf_dep(ji,jj) = 0.
-            DO jk = 1, rnf_mod_dep(ji,jj)                          ! recalculates rnf_dep to be the depth 
-               rnf_dep(ji,jj) = rnf_dep(ji,jj) + fse3t(ji,jj,jk)    ! in metres to the bottom of the relevant grid box 
-            ENDDO
-         ENDDO
-      ENDDO
-      ! 
 
       !                                   ! ========================

branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/TRA/trabbc.F90

@@ -36 +36 @@
    REAL(wp) ::   rn_geoflx_cst = 86.4e-3      ! Constant value of geothermal heat flux
 
-   INTEGER , DIMENSION(jpi,jpj) ::   nbotlevt   ! ocean bottom level index at T-pt
-   REAL(wp), DIMENSION(jpi,jpj) ::   qgh_trd0   ! geothermal heating trend
+   INTEGER , DIMENSION(jpi,jpj)         ::   nbotlevt   ! ocean bottom level index at T-pt
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qgh_trd0   ! geothermal heating trend
  
    !! * Substitutions

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

@@ -138 +138 @@
 
 
-   SUBROUTINE tra_bbl_dif( ptrab, ptraa, kjpt )
+   SUBROUTINE tra_bbl_dif( ptb, pta, kjpt )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_bbl_dif  ***
@@ -155 +155 @@
       !!      convection is satified)
       !!
-      !! ** Action  :   ptraa   increased by the bbl diffusive trend
+      !! ** Action  :   pta   increased by the bbl diffusive trend
       !!
       !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
@@ -161 +161 @@
       !!----------------------------------------------------------------------  
       INTEGER                              , INTENT(in   ) ::   kjpt    ! number of tracers
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptrab   ! before and now tracer fields
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   ptraa   ! tracer trend 
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb   ! before and now tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta   ! tracer trend 
       !!
       INTEGER  ::   ji, jj, jn   ! dummy loop indices
       INTEGER  ::   ik           ! local integers
       REAL(wp) ::   zbtr         ! local scalars
+      REAL(wp), DIMENSION(jpi,jpj) ::   zptb   ! tracer trend 
       !!----------------------------------------------------------------------
       !
-      !                                                   ! ===========
       DO jn = 1, kjpt                                     ! tracer loop
          !                                                ! ===========
+#  if defined key_vectopt_loop
+         DO jj = 1, 1   ! vector opt. (forced unrolling)
+            DO ji = 1, jpij
+#else
+         DO jj = 1, jpj
+            DO ji = 1, jpi 
+#endif
+               ik = mbkt(ji,jj)                        ! bottom T-level index
+               zptb(ji,jj) = ptb(ji,jj,ik,jn)              ! bottom before T and S
+            END DO
+         END DO
+         !                                                ! Compute the trend
 #  if defined key_vectopt_loop
          DO jj = 1, 1   ! vector opt. (forced unrolling)
@@ -181 +193 @@
                ik = mbkt(ji,jj)                            ! bottom T-level index
                zbtr = e1e2t_r(ji,jj)  / fse3t(ji,jj,ik)
-               ptraa(ji,jj,ik,jn) = ptraa(ji,jj,ik,jn)                                                         &
-                  &               + (   ahu_bbl(ji,jj) * ( ptrab(ji+1,jj  ,ik,jn) - ptrab(ji  ,jj  ,ik,jn) )   &
-                  &                   - ahu_bbl(ji,jj) * ( ptrab(ji  ,jj  ,ik,jn) - ptrab(ji-1,jj  ,ik,jn) )   &
-                  &                   + ahv_bbl(ji,jj) * ( ptrab(ji  ,jj+1,ik,jn) - ptrab(ji  ,jj  ,ik,jn) )   &
-                  &                   - ahv_bbl(ji,jj) * ( ptrab(ji  ,jj  ,ik,jn) - ptrab(ji  ,jj-1,ik,jn) )   ) * zbtr
+               pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                         &
+                  &               + (   ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )   &
+                  &                   - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )   &
+                  &                   + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )   &
+                  &                   - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )   ) * zbtr
             END DO
          END DO
@@ -194 +206 @@
    
 
-   SUBROUTINE tra_bbl_adv( ptrab, ptraa, kjpt )
+   SUBROUTINE tra_bbl_adv( ptb, pta, kjpt )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE trc_bbl  ***
@@ -212 +224 @@
       !!----------------------------------------------------------------------  
       INTEGER                              , INTENT(in   ) ::   kjpt    ! number of tracers
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptrab   ! before and now tracer fields
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   ptraa   ! tracer trend 
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb   ! before and now tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta   ! tracer trend 
       !!
       INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
@@ -240 +252 @@
                   !                                               ! up  -slope T-point (shelf bottom point)
                   zbtr = e1e2t_r(iis,jj) / fse3t(iis,jj,ikus)
-                  ztra = zu_bbl * ( ptrab(iid,jj,ikus,jn) - ptrab(iis,jj,ikus,jn) ) * zbtr
-                  ptraa(iis,jj,ikus,jn) = ptraa(iis,jj,ikus,jn) + ztra
+                  ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr
+                  pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra
                   !                   
                   DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
                      zbtr = e1e2t_r(iid,jj) / fse3t(iid,jj,jk)
-                     ztra = zu_bbl * ( ptrab(iid,jj,jk+1,jn) - ptrab(iid,jj,jk,jn) ) * zbtr
-                     ptraa(iid,jj,jk,jn) = ptraa(iid,jj,jk,jn) + ztra
+                     ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr
+                     pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra
                   END DO
                   ! 
                   zbtr = e1e2t_r(iid,jj) / fse3t(iid,jj,ikud)
-                  ztra = zu_bbl * ( ptrab(iis,jj,ikus,jn) - ptrab(iid,jj,ikud,jn) ) * zbtr
-                  ptraa(iid,jj,ikud,jn) = ptraa(iid,jj,ikud,jn) + ztra
+                  ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr
+                  pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra
                ENDIF
                !
@@ -262 +274 @@
                   ! up  -slope T-point (shelf bottom point)
                   zbtr = e1e2t_r(ji,ijs) / fse3t(ji,ijs,ikvs)
-                  ztra = zv_bbl * ( ptrab(ji,ijd,ikvs,jn) - ptrab(ji,ijs,ikvs,jn) ) * zbtr
-                  ptraa(ji,ijs,ikvs,jn) = ptraa(ji,ijs,ikvs,jn) + ztra
+                  ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr
+                  pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra
                   !                   
                   DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
                      zbtr = e1e2t_r(ji,ijd) / fse3t(ji,ijd,jk)
-                     ztra = zv_bbl * ( ptrab(ji,ijd,jk+1,jn) - ptrab(ji,ijd,jk,jn) ) * zbtr
-                     ptraa(ji,ijd,jk,jn) = ptraa(ji,ijd,jk,jn)  + ztra
+                     ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr
+                     pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn)  + ztra
                   END DO
                   !                                               ! down-slope T-point (deep bottom point)
                   zbtr = e1e2t_r(ji,ijd) / fse3t(ji,ijd,ikvd)
-                  ztra = zv_bbl * ( ptrab(ji,ijs,ikvs,jn) - ptrab(ji,ijd,ikvd,jn) ) * zbtr
-                  ptraa(ji,ijd,ikvd,jn) = ptraa(ji,ijd,ikvd,jn) + ztra
+                  ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr
+                  pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra
                ENDIF
             END DO
@@ -370 +382 @@
 #endif
             ik = mbkt(ji,jj)                        ! bottom T-level index
-            ztb (ji,jj) = tsb(ji,jj,ik,jp_tem)      ! bottom before T and S
-            zsb (ji,jj) = tsb(ji,jj,ik,jp_sal)
+            ztb (ji,jj) = tsb(ji,jj,ik,jp_tem) * tmask(ji,jj,1)      ! bottom before T and S
+            zsb (ji,jj) = tsb(ji,jj,ik,jp_sal) * tmask(ji,jj,1)
             zdep(ji,jj) = fsdept_0(ji,jj,ik)        ! bottom T-level reference depth
             !
@@ -575 +587 @@
       ahv_bbl_0(:,:) = rn_ahtbbl * e1v(:,:) * e3v_bbl_0(:,:) / e2v(:,:)  * vmask(:,:,1)
 
+
       IF( cp_cfg == "orca" ) THEN   !* ORCA configuration : regional enhancement of ah_bbl
          !

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

@@ -15 +15 @@
    !!            3.0  !  2008-06  (G. Madec)  time stepping always done in trazdf
    !!            3.1  !  2009-02  (G. Madec, R. Benshila)  re-introduce the vvl option
-   !!            3.3  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  semi-implicit hpg with asselin filter + modified LF-RA
+   !!             -   !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
-   !!   tra_nxt       : time stepping on temperature and salinity
-   !!   tra_nxt_fix   : time stepping on temperature and salinity : fixed    volume case
-   !!   tra_nxt_vvl   : time stepping on temperature and salinity : variable volume case
+   !!   tra_nxt       : time stepping on tracers
+   !!   tra_nxt_fix   : time stepping on tracers : fixed    volume case
+   !!   tra_nxt_vvl   : time stepping on tracers : variable volume case
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers variables
    USE dom_oce         ! ocean space and time domain variables 
+   USE sbc_oce         ! surface boundary condition: ocean
    USE zdf_oce         ! ???
    USE domvvl          ! variable volume
@@ -38 +40 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE prtctl          ! Print control
+   USE traqsr          ! penetrative solar radiation (needed for nksr)
    USE traswp          ! swap array
    USE agrif_opa_update
@@ -49 +52 @@
    PUBLIC   tra_nxt_vvl   ! to be used in trcnxt
 
+   REAL(wp)                 ::   rbcp            ! Brown & Campana parameters for semi-implicit hpg
    REAL(wp), DIMENSION(jpk) ::   r2dt   ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler)
 
@@ -96 +100 @@
          IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap'
          IF(lwp) WRITE(numout,*) '~~~~~~~'
+         !
+         rbcp    = 0.25 * (1. + atfp) * (1. + atfp) * ( 1. - atfp)       ! Brown & Campana parameter for semi-implicit hpg
       ENDIF
 
@@ -107 +113 @@
 #endif
 
-#if defined key_obc
+#if defined key_obc 
       IF( lk_obc )   CALL obc_tra( kt )  ! OBC open boundaries
 #endif
@@ -114 +120 @@
 #endif
 #if defined key_agrif
-      CALL Agrif_tra                     ! AGRIF zoom boundaries
+      CALL Agrif_tra                   ! AGRIF zoom boundaries
 #endif
 
@@ -133 +139 @@
 
       ! Leap-Frog + Asselin filter time stepping
-      IF( lk_vvl )   THEN   ;   CALL tra_nxt_vvl( kt, tsb, tsn, tsa, jpts )  ! variable volume level (vvl)     
-      ELSE                  ;   CALL tra_nxt_fix( kt, tsb, tsn, tsa, jpts )  ! fixed    volume level 
+      IF( lk_vvl )   THEN   ;   CALL tra_nxt_vvl( kt, 'TRA', tsb, tsn, tsa, jpts )  ! variable volume level (vvl)     
+      ELSE                  ;   CALL tra_nxt_fix( kt,        tsb, tsn, tsa, jpts )  ! fixed    volume level 
       ENDIF
 
@@ -176 +182 @@
       !!              - swap tracer fields to prepare the next time_step.
       !!                This can be summurized for tempearture as:
-      !!             ztm = (ta+2tn+tb)/4       ln_dynhpg_imp = T
-      !!             ztm = 0                   otherwise
+      !!             ztm = tn + rbcp * [ta -2 tn + tb ]       ln_dynhpg_imp = T
+      !!             ztm = 0                                   otherwise
+      !!                   with rbcp=1/4 * (1-atfp^4) / (1-atfp)
       !!             tb  = tn + atfp*[ tb - 2 tn + ta ]
-      !!             tn  = ta 
+      !!             tn  = ta  
       !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call)
       !!
@@ -185 +192 @@
       !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
       !!----------------------------------------------------------------------
-      INTEGER , INTENT(in   )                               ::  kt            ! ocean time-step index
-      INTEGER , INTENT(in   )                               ::  kjpt            ! number of tracers
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb  ! before tracer fields
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn  ! now tracer fields
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta         ! tracer trend
+      INTEGER , INTENT(in   )                               ::  kt       ! ocean time-step index
+      INTEGER , INTENT(in   )                               ::  kjpt     ! number of tracers
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb      ! before tracer fields
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn      ! now tracer fields
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta      ! tracer trend
       !!
       INTEGER  :: ji, jj, jk, jn   ! dummy loop indices
-      REAL(wp) :: ztm, ztf     ! temporary scalars
+      REAL(wp) :: ztd, ztm         ! temporary scalars
       !!----------------------------------------------------------------------
 
@@ -205 +212 @@
          !                                             ! (only swap)
          DO jn = 1, kjpt
-            DO jk = 1, jpkm1                              
+            DO jk = 1, jpkm1
                ptn(:,:,jk,jn) = pta(:,:,jk,jn)     ! ptb <-- ptn
             END DO
@@ -219 +226 @@
                   DO jj = 1, jpj
                      DO ji = 1, jpi
-                        ztm = 0.25 * ( pta(ji,jj,jk,jn) + 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) )  ! mean pt
-                        ztf = atfp * ( pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) )  ! Asselin filter on pt 
-                        ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + ztf                                      ! ptb <-- filtered ptn 
-                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                                            ! ptn <-- pta
-                        pta(ji,jj,jk,jn) = ztm                                                           ! pta <-- mean pt
+                        ztd = pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn)    !  time laplacian on tracers
+                        ztm = ptn(ji,jj,jk,jn) + rbcp * ztd                                 !  used for both Asselin and Brown & Campana filters
+                        !
+                        ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + atfp * ztd                    ! ptb <-- filtered ptn 
+                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                                 ! ptn <-- pta
+                        pta(ji,jj,jk,jn) = ztm                                              ! pta <-- Brown & Campana average
                      END DO
                   END DO
@@ -235 +243 @@
                   DO jj = 1, jpj
                      DO ji = 1, jpi
-                        ztf = atfp * ( pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) )       ! Asselin filter on t 
-                        ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + ztf                                     ! ptb <-- filtered ptn 
-                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                                           ! ptn <-- pta
+                        ztd = pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn)    ! time laplacian on tracers
+                        ! 
+                        ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + atfp * ztd                    ! ptb <-- filtered ptn 
+                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                                 ! ptn <-- pta
                      END DO
                   END DO
                END DO
             END DO
-            !
          ENDIF
          !
@@ -250 +258 @@
 
 
-   SUBROUTINE tra_nxt_vvl( kt, ptb, ptn, pta, kjpt )
+   SUBROUTINE tra_nxt_vvl( kt, cdtype, ptb, ptn, pta, kjpt )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE tra_nxt_vvl  ***
@@ -263 +271 @@
       !!              - swap tracer fields to prepare the next time_step.
       !!                This can be summurized for tempearture as:
-      !!             ztm = (e3t_a*ta+2*e3t_n*tn+e3t_b*tb)   ln_dynhpg_imp = T
-      !!                  /(e3t_a   +2*e3t_n   +e3t_b   )    
-      !!             ztm = 0                                otherwise
+      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
+      !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )   
+      !!             ztm = 0                                                       otherwise
       !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
       !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
@@ -274 +282 @@
       !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
       !!----------------------------------------------------------------------
-      INTEGER , INTENT(in   )                               ::  kt            ! ocean time-step index
-      INTEGER , INTENT(in   )                               ::  kjpt            ! number of tracers
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb  ! before tracer fields
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn  ! now tracer fields
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta         ! tracer trend
+      INTEGER         , INTENT(in   )                               ::  kt       ! ocean time-step index
+      CHARACTER(len=3), INTENT(in   )                               ::  cdtype   ! =TRA or TRC (tracer indicator)
+      INTEGER         , INTENT(in   )                               ::  kjpt     ! number of tracers
+      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb      ! before tracer fields
+      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn      ! now tracer fields
+      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta      ! tracer trend
       !!     
-      INTEGER  ::   ji, jj, jk, jn             ! dummy loop indices
-      REAL(wp) ::   ztm , ztc_f , ztf , ztca, ztcn, ztcb   ! temporary scalar
-      REAL(wp) ::   ze3mr, ze3fr                           !    -         -
-      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f         !    -         -
+      INTEGER  ::   ji, jj, jk, jn                 ! dummy loop indices
+      REAL(wp) ::   ztc_a , ztc_n , ztc_b          ! temporary scalar
+      REAL(wp) ::   ztc_f , ztc_d , ztc_m          !    -         -
+      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a         !    -         -
+      REAL(wp) ::   ze3t_f, ze3t_d, ze3t_m         !    -         -
+      REAL     ::   zfact1, zfact2                 !    -         -
       !!----------------------------------------------------------------------
 
@@ -293 +304 @@
       !
       !
-      IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step
-         !                                             ! (only swap)
+      IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step
+         !                                           ! (only swap)
          DO jn = 1, kjpt
-            DO jk = 1, jpkm1                              
-               ptn(:,:,jk,jn) = pta(:,:,jk,jn)         ! ptb <-- ptn
+            DO jk = 1, jpkm1
+               ptn(:,:,jk,jn) = pta(:,:,jk,jn)       ! ptb <-- ptn
             END DO
          END DO
          !                                              
-      ELSE                                              ! general case (Leapfrog + Asselin filter
+      ELSE                                           ! general case (Leapfrog + Asselin filter)
          !
-         !                                              ! ----------------------- !
-         IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg case !
-            !                                           ! ----------------------- !
-            DO jn = 1, kjpt                               
+         !                                           ! ----------------------- !
+         IF( ln_dynhpg_imp ) THEN                    ! semi-implicite hpg case !
+            !                                        ! ----------------------- !
+            DO jn = 1, kjpt                          
                DO jk = 1, jpkm1
                   DO jj = 1, jpj
@@ -314 +325 @@
                         ze3t_a = fse3t_a(ji,jj,jk)
                         !                                         ! tracer content at Before, now and after
-                        ztcb = ptb(ji,jj,jk,jn) *  ze3t_b  
-                        ztcn = ptn(ji,jj,jk,jn) *  ze3t_n  
-                        ztca = pta(ji,jj,jk,jn) *  ze3t_a  
-                        !
-                        !                                         ! Asselin filter on thickness and tracer content
-                        ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b )
-                        ztc_f  = atfp * ( ztca   - 2.* ztcn   + ztcb   ) 
-                        !
-                        !                                         ! filtered tracer including the correction 
-                        ze3fr = 1.e0  / ( ze3t_n + ze3t_f )
-                        ztf   = ze3fr * ( ztcn   + ztc_f  )
-                        !                                         ! mean thickness and tracer
-                        ze3mr = 1.e0  / ( ze3t_a + 2.* ze3t_n + ze3t_b )
-                        ztm   = ze3mr * ( ztca   + 2.* ztcn   + ztcb   )
-                        !!gm mean e3t have to be saved and used in dynhpg  or it can be recomputed in dynhpg !!
-                        !!gm                 e3t_m(ji,jj,jk) = 0.25 / ze3mr
+                        ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b  
+                        ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n  
+                        ztc_a  = pta(ji,jj,jk,jn) * ze3t_a  
+                        !                                         ! Time laplacian on tracer contents
+                        !                                         ! used for both Asselin and Brown & Campana filters
+                        ze3t_d =  ze3t_b - 2.* ze3t_n + ze3t_a 
+                        ztc_d  =  ztc_b  - 2.* ztc_n  + ztc_a  
+                        !                                         ! Asselin Filter on thicknesses and tracer contents
+                        ztc_f  = ztc_n  + atfp * ztc_d
+                        ztc_m  = ztc_n  + rbcp * ztc_d
+                        !                                         
+                        ze3t_f = 1.0 / ( ze3t_n + atfp * ze3t_d )
+                        ze3t_m = 1.0 / ( ze3t_n + rbcp * ze3t_d )
                         !                                         ! swap of arrays
-                        ptb(ji,jj,jk,jn) = ztf                             ! ptb <-- ptn + filter
-                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)              ! ptn <-- pta
-                        pta(ji,jj,jk,jn) = ztm                             ! pta <-- mean t
+                        ptb(ji,jj,jk,jn) = ztc_f * ze3t_f              ! ptb <-- ptn filtered
+                        pta(ji,jj,jk,jn) = ztc_m * ze3t_m              ! pta <-- Brown & Campana average
+                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)            ! ptn <-- pta
                      END DO
                   END DO
                END DO
             END DO
-            !                                        ! ----------------------- !
-         ELSE                                        !    explicit hpg case    !
-            !                                        ! ----------------------- !
-            DO jn = 1, kjpt      
-               DO jk = 1, jpkm1
-                  DO jj = 1, jpj
-                     DO ji = 1, jpi
-                        ze3t_b = fse3t_b(ji,jj,jk)
-                        ze3t_n = fse3t_n(ji,jj,jk)
-                        ze3t_a = fse3t_a(ji,jj,jk)
-                        !                                         ! tracer content at Before, now and after
-                        ztcb = ptb(ji,jj,jk,jn) *  ze3t_b  
-                        ztcn = ptn(ji,jj,jk,jn) *  ze3t_n   
-                        ztca = pta(ji,jj,jk,jn) *  ze3t_a  
-                        !
-                        !                                         ! Asselin filter on thickness and tracer content
-                        ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b )
-                        ztc_f  = atfp * ( ztca   - 2.* ztcn   + ztcb   ) 
-                        !
-                        !                                         ! filtered tracer including the correction 
-                        ze3fr = 1.e0 / ( ze3t_n + ze3t_f )
-                        ztf   =  ( ztcn  + ztc_f ) * ze3fr
-                        !                                         ! swap of arrays
-                        ptb(ji,jj,jk,jn) = ztf                    ! tb <-- tn filtered
-                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)       ! tn <-- ta
+            !                                           ! ----------------------- !
+         ELSE                                           !    explicit hpg case    !
+            !                                           ! ----------------------- !
+            IF( cdtype == 'TRA' ) THEN
+               !              
+               DO jn = 1, kjpt      
+                  DO jk = 1, jpkm1
+                     zfact1 = atfp * rdttra(jk)
+                     zfact2 = zfact1 / rau0
+                     DO jj = 1, jpj
+                        DO ji = 1, jpi
+                           ze3t_b = fse3t_b(ji,jj,jk)
+                           ze3t_n = fse3t_n(ji,jj,jk)
+                           ze3t_a = fse3t_a(ji,jj,jk)
+                           !                                         ! tracer content at Before, now and after
+                           ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b
+                           ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n
+                           ztc_a  = pta(ji,jj,jk,jn) * ze3t_a
+                           !
+                           ze3t_f = ze3t_n + atfp * ( ze3t_b - 2. * ze3t_n + ze3t_a )
+                           ztc_f  = ztc_n  + atfp * ( ztc_a  - 2. * ztc_n  + ztc_b )
+
+                           IF( jk == 1 ) THEN
+                               ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) )
+                               IF( jn == jp_tem )   ztc_f  = ztc_f  - zfact1 * ( sbc_hc_n(ji,jj) - sbc_hc_b(ji,jj) )
+                               IF( jn == jp_sal )   ztc_f  = ztc_f  - zfact1 * ( sbc_sc_n(ji,jj) - sbc_sc_b(ji,jj) )
+                           ENDIF
+                           IF( jn == jp_tem .AND. ln_traqsr .AND. jk <= nksr )  &
+                           &                        ztc_f  = ztc_f  - zfact1 * ( qsr_hc_n(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
+
+                           ze3t_f = 1.e0 / ze3t_f
+                           ptb(ji,jj,jk,jn) = ztc_f * ze3t_f                           ! tb <-- tn filtered
+                           ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                             ! tn <-- ta
+                        END DO
                      END DO
                   END DO
                END DO
-            END DO
+               !
+            ELSE IF( cdtype == 'TRC' ) THEN
+               !              
+               DO jn = 1, kjpt      
+                  DO jk = 1, jpkm1
+                     DO jj = 1, jpj
+                        DO ji = 1, jpi
+                           ze3t_b = fse3t_b(ji,jj,jk)
+                           ze3t_n = fse3t_n(ji,jj,jk)
+                           ze3t_a = fse3t_a(ji,jj,jk)
+                           !                                         ! tracer content at Before, now and after
+                           ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b
+                           ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n
+                           ztc_a  = pta(ji,jj,jk,jn) * ze3t_a
+                           !
+                           ze3t_f = ze3t_n + atfp * ( ze3t_b - 2. * ze3t_n + ze3t_a )
+                           ztc_f  = ztc_n  + atfp * ( ztc_a  - 2. * ztc_n  + ztc_b  )
+
+                           ze3t_f = 1.e0 / ze3t_f
+                           ptb(ji,jj,jk,jn) = ztc_f * ze3t_f                           ! tb <-- tn filtered
+                           ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                             ! tn <-- ta
+                        END DO
+                     END DO
+                  END DO
+               END DO
+               !
+            END IF
             !
          ENDIF

branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/TRA/traqsr.F90

@@ -27 +27 @@
    USE iom             ! I/O manager
    USE fldread         ! read input fields
+   USE restart         ! ocean restart
 
    IMPLICIT NONE
@@ -47 +48 @@
    ! Module variables
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read)
-   INTEGER ::   nksr   ! levels below which the light cannot penetrate ( depth larger than 391 m)
+   INTEGER, PUBLIC ::   nksr              ! levels below which the light cannot penetrate ( depth larger than 391 m)
    REAL(wp), DIMENSION(3,61) ::   rkrgb   !: tabulated attenuation coefficients for RGB absorption
 
@@ -95 +96 @@
       REAL(wp) ::   zchl, zcoef, zsi0r   ! temporary scalars
       REAL(wp) ::   zc0, zc1, zc2, zc3   !    -         -
+      REAL(wp) ::   z1_e3t, zfact        !    -         -
       REAL(wp), DIMENSION(jpi,jpj)     ::   zekb, zekg, zekr            ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze0, ze1 , ze2, ze3, zea    ! 3D workspace
@@ -111 +113 @@
          ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = 0.
       ENDIF
+
+      !                                        Set before qsr tracer content field
+      !                                        ***********************************
+      IF( kt == nit000 ) THEN                     ! Set the forcing field at nit000 - 1
+         !                                        ! -----------------------------------
+         IF( ln_rstart .AND.    &                    ! Restart: read in restart file
+              & iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN
+            IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field red in the restart file'
+            zfact = 0.5e0
+            CALL iom_get( numror, jpdom_autoglo, 'qsr_hc_b', qsr_hc_b )   ! before heat content trend due to Qsr flux
+         ELSE                                           ! No restart or restart not found: Euler forward time stepping
+            zfact = 1.e0
+            qsr_hc_b(:,:,:) = 0.e0
+         ENDIF
+      ELSE                                        ! Swap of forcing field
+         !                                        ! ---------------------
+         zfact = 0.5e0
+         qsr_hc_b(:,:,:) = qsr_hc_n(:,:,:)
+      ENDIF
+      !                                        Compute now qsr tracer content field
+      !                                        ************************************
       
       !                                           ! ============================================== !
@@ -118 +141 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ta(ji,jj,jk) = ta(ji,jj,jk) + ro0cpr * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) / fse3t(ji,jj,jk) 
+                  qsr_hc_n(ji,jj,jk) =  ro0cpr * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) / fse3t(ji,jj,jk) 
                END DO
             END DO
@@ -175 +198 @@
                !
                DO jk = 1, nksr                                        ! compute and add qsr trend to ta
-                  tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t(:,:,jk)
+                  qsr_hc_n(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) )
                END DO
                zea(:,:,nksr+1:jpk) = 0.e0     ! below 400m set to zero
@@ -182 +205 @@
             ELSE                                                 !*  Constant Chlorophyll
                DO jk = 1, nksr
-                  tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + etot3(:,:,jk) * qsr(:,:)
+                  qsr_hc_n(:,:,jk) =  etot3(:,:,jk) * qsr(:,:)
                END DO
             ENDIF
@@ -194 +217 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + etot3(ji,jj,jk) * qsr(ji,jj)
+                     qsr_hc_n(ji,jj,jk) =  etot3(ji,jj,jk) * qsr(ji,jj)
                   END DO
                END DO
@@ -200 +223 @@
             !
          ENDIF
+         !
+      ENDIF
+      !                                        Add to the general trend
+      !                                        ************************
+      DO jk = 1, nksr
+         DO jj = 2, jpjm1 
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               z1_e3t = zfact / fse3t(ji,jj,jk)
+               tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + ( qsr_hc_b(ji,jj,jk) + qsr_hc_n(ji,jj,jk) ) * z1_e3t
+            END DO
+         END DO
+      END DO
+      !
+      IF( lrst_oce ) THEN   !                  Write in the ocean restart file
+         !                                     *******************************
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'qsr tracer content forcing field written in ocean restart file ',   &
+            &                    'at it= ', kt,' date= ', ndastp
+         IF(lwp) WRITE(numout,*) '~~~~'
+         CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b', qsr_hc_n )
          !
       ENDIF

branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/TRA/trasbc.F90

@@ -7 +7 @@
    !!            8.2  !  2001-02  (D. Ludicone)  sea ice and free surface
    !!  NEMO      1.0  !  2002-06  (G. Madec)  F90: Free form and module
-   !!            3.3  !  2010-09  (C. Ethe, G. Madec) Merge TRA-TRC
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
+   !!             -   !  2010-09  (C. Ethe, G. Madec) Merge TRA-TRC
    !!----------------------------------------------------------------------
 
@@ -22 +23 @@
    USE in_out_manager  ! I/O manager
    USE prtctl          ! Print control
+   USE restart         ! ocean restart
    USE sbcrnf          ! River runoff  
    USE sbcmod          ! ln_rnf  
+   USE iom
 
    IMPLICIT NONE
@@ -103 +106 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !!
-      INTEGER  ::   ji, jj, jk      ! dummy loop indices  
-      REAL(wp) ::   zta, zsa        ! local scalars, adjustment to temperature and salinity  
-      REAL(wp) ::   zata, zasa      ! local scalars, calculations of automatic change to temp & sal due to vvl (done elsewhere)  
+      INTEGER  ::   ji, jj, jk           ! dummy loop indices  
+      REAL(wp) ::   zta, zsa, zrnf       ! local scalars, adjustment to temperature and salinity  
       REAL(wp) ::   zsrau, zse3t, zdep   ! local scalars, 1/density, 1/height of box, 1/height of effected water column  
       REAL(wp) ::   zdheat, zdsalt       ! total change of temperature and salinity  
+      REAL(wp) ::   zfact, z1_e3t        ! 
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ztrdt, ztrds
       !!----------------------------------------------------------------------
@@ -118 +121 @@
 
       zsrau = 1. / rau0             ! initialization
-#if defined key_zco
-      zse3t = 1. / e3t_0(1)
-#endif
 
       IF( l_trdtra )   THEN                    !* Save ta and sa trends
@@ -127 +127 @@
       ENDIF
 
-      IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration
-
-      ! Concentration dilution effect on (t,s) due to evapouration, precipitation and qns, but not river runoff  
+!!gm      IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration
+      IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration
+         qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns
+         qsr(:,:) = 0.e0                     ! qsr set to zero
+      ENDIF
+
+      !                                          Set before sbc tracer content fields
+      !                                          ************************************
+      IF( kt == nit000 ) THEN                      ! Set the forcing field at nit000 - 1
+         !                                         ! -----------------------------------
+         IF( ln_rstart .AND.    &                     ! Restart: read in restart file
+              & iom_varid( numror, 'sbc_hc_b', ldstop = .FALSE. ) > 0 ) THEN
+            IF(lwp) WRITE(numout,*) '          nit000-1 surface tracer content forcing fields red in the restart file'
+            zfact = 0.5e0
+            CALL iom_get( numror, jpdom_autoglo, 'sbc_hc_b', sbc_hc_b )   ! before heat content sbc trend
+            CALL iom_get( numror, jpdom_autoglo, 'sbc_sc_b', sbc_sc_b )   ! before salt content sbc trend
+         ELSE                                         ! No restart or restart not found: Euler forward time stepping
+            zfact = 1.e0
+            sbc_hc_b(:,:) = 0.e0
+            sbc_sc_b(:,:) = 0.e0
+         ENDIF
+      ELSE                                         ! Swap of forcing fields
+         !                                         ! ----------------------
+         zfact = 0.5e0
+         sbc_hc_b(:,:) = sbc_hc_n(:,:)
+         sbc_sc_b(:,:) = sbc_sc_n(:,:)
+      ENDIF
+      !                                          Compute now sbc tracer content fields
+      !                                          *************************************
+
+                                                   ! Concentration dilution effect on (t,s) due to  
+                                                   ! evaporation, precipitation and qns, but not river runoff 
+                                               
+      IF( lk_vvl ) THEN                            ! Variable Volume case
+         DO jj = 2, jpj
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ! temperature : heat flux + cooling/heating effet of EMP flux
+               sbc_hc_n(ji,jj) = ro0cpr * qns(ji,jj) - zsrau * emp(ji,jj) * tsn(ji,jj,1,jp_tem)
+               ! concent./dilut. effect due to sea-ice melt/formation and (possibly) SSS restoration
+               sbc_sc_n(ji,jj) = ( emps(ji,jj) - emp(ji,jj) ) * zsrau * tsn(ji,jj,1,jp_sal)
+            END DO
+         END DO
+      ELSE                                         ! Constant Volume case
+         DO jj = 2, jpj
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ! temperature : heat flux
+               sbc_hc_n(ji,jj) = ro0cpr * qns(ji,jj)
+               ! salinity    : salt flux + concent./dilut. effect (both in emps)
+               sbc_sc_n(ji,jj) = zsrau * emps(ji,jj) * tsn(ji,jj,1,jp_sal)
+            END DO
+         END DO
+      ENDIF
+                                                   ! Concentration dilution effect on (t,s) due to  
+                                                   ! river runoff without T, S and depth attributes
+      IF( ln_rnf ) THEN 
+         !
+         IF( lk_vvl ) THEN                            ! Variable Volume case 
+            !                                         ! cooling/heating effect of runoff & No salinity concent./dilut. effect
+            DO jj = 2, jpj
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  sbc_hc_n(ji,jj) = sbc_hc_n(ji,jj) + zsrau * rnf(ji,jj) * tsn(ji,jj,1,jp_tem)
+               END DO
+            END DO
+         ELSE                                         ! Constant Volume case
+            !                                         ! concent./dilut. effect only
+            DO jj = 2, jpj 
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  sbc_sc_n(ji,jj) = sbc_sc_n(ji,jj) - zsrau * rnf(ji,jj) * tsn(ji,jj,1,jp_sal)
+               END DO
+            END DO
+         ENDIF
+         !  
+      ENDIF
+      !                                          Add to the general trend
+      !                                          ************************
       DO jj = 2, jpj
          DO ji = fs_2, fs_jpim1   ! vector opt.
-#if ! defined key_zco
-            zse3t = 1. / fse3t(ji,jj,1)
-#endif
-            IF( lk_vvl) THEN
-               ! temperature : heat flux and heat content of EMP flux
-               zta = ( ro0cpr * qns(ji,jj) - emp(ji,jj) * zsrau * tsn(ji,jj,1,jp_tem) ) * zse3t
-               ! Salinity    : concent./dilut. effect due to sea-ice melt/formation and (possibly) SSS restoration
-               zsa = ( emps(ji,jj) - emp(ji,jj) ) * zsrau * tsn(ji,jj,1,jp_sal) * zse3t
-            ELSE
-               zta =  ro0cpr * qns(ji,jj) * zse3t                         ! temperature : heat flux 
-               zsa =  emps(ji,jj) * zsrau * tsn(ji,jj,1,jp_sal) * zse3t   ! salinity :  concent./dilut. effect 
-            ENDIF
-            tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + zta                  ! add the trend to the general tracer trend
-            tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + zsa
+            z1_e3t = zfact / fse3t(ji,jj,1)
+            tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + ( sbc_hc_b(ji,jj) + sbc_hc_n(ji,jj) ) * z1_e3t
+            tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + ( sbc_sc_b(ji,jj) + sbc_sc_n(ji,jj) ) * z1_e3t
          END DO
       END DO
+      !                                          Write in the ocean restart file
+      !                                          *******************************
+      IF( lrst_oce ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'sbc : ocean surface tracer content forcing fields written in ocean restart file ',   &
+            &                    'at it= ', kt,' date= ', ndastp
+         IF(lwp) WRITE(numout,*) '~~~~'
+         CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_hc_n )
+         CALL iom_rstput( kt, nitrst, numrow, 'sbc_sc_b', sbc_sc_n )
+      ENDIF
 
       IF( ln_rnf .AND. ln_rnf_att ) THEN        ! Concentration / dilution effect on (t,s) due to river runoff  
+        !
         DO jj = 1, jpj  
            DO ji = 1, jpi  
-              zdep = 1. / rnf_dep(ji,jj)  
-              zse3t= 1. / fse3t(ji,jj,1)  
-              IF( rnf_tem(ji,jj) == -999 )   rnf_tem(ji,jj) = tsn(ji,jj,1,jp_tem)   ! if not specified set runoff temp to be sst  
+              zdep  = 1. / rnf_dep(ji,jj)  
+              zse3t = 1. / fse3t(ji,jj,1)  
+              rnf_dep(ji,jj) = 0.e0
+              DO jk = 1, rnf_mod_dep(ji,jj)                          ! recalculates rnf_dep to be the depth  
+                rnf_dep(ji,jj) = rnf_dep(ji,jj) + fse3t(ji,jj,jk)    ! in metres to the bottom of the relevant grid box  
+              END DO
+              IF( rnf_tmp(ji,jj) == -999 )   rnf_tmp(ji,jj) = tsn(ji,jj,1,jp_tem)   ! if not specified set runoff temp to be sst  
   
               IF( rnf(ji,jj) > 0.e0 ) THEN  
-  
+                !
+                zrnf = rnf(ji,jj) * zsrau * zdep
                 IF( lk_vvl ) THEN  
                   ! indirect flux, concentration or dilution effect : force a dilution effect in all levels 
@@ -163 +240 @@
                   zdsalt = 0.e0  
                   DO jk = 1, rnf_mod_dep(ji,jj)  
-                    zta = -tsn(ji,jj,jk,jp_tem) * rnf(ji,jj) * zsrau * zdep  
-                    zsa = -tsn(ji,jj,jk,jp_sal) * rnf(ji,jj) * zsrau * zdep  
+                    zta = -tsn(ji,jj,jk,jp_tem) * zrnf  
+                    zsa = -tsn(ji,jj,jk,jp_sal) * zrnf  
                     tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta        ! add the trend to the general tracer trend
                     tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa
@@ -171 +248 @@
                   END DO  
                   ! negate this total change in heat and salt content from top level    !!gm I don't understand this
-                  zta = -zdheat * zse3t  
-                  zsa = -zdsalt * zse3t  
-                  tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + zta            ! add the trend to the general tracer trend
-                  tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + zsa
-    
-                  ! direct flux  
-                  zta = rnf_tem(ji,jj) * rnf(ji,jj) * zsrau * zdep  
-                  zsa = rnf_sal(ji,jj) * rnf(ji,jj) * zsrau * zdep  
+                  tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) - zdheat * zse3t
+                  tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) - zdsalt * zse3t
     
                   DO jk = 1, rnf_mod_dep(ji,jj)  
-                    tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta        ! add the trend to the general tracer trend
-                    tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa
+                    tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + rnf_tmp(ji,jj) * zrnf
+                    tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + rnf_sal(ji,jj) * zrnf
                   END DO  
                 ELSE  
                   DO jk = 1, rnf_mod_dep(ji,jj)  
-                    zta = ( rnf_tem(ji,jj) - tsn(ji,jj,jk,jp_tem) ) * rnf(ji,jj) * zsrau * zdep  
-                    zsa = ( rnf_sal(ji,jj) - tsn(ji,jj,jk,jp_sal) ) * rnf(ji,jj) * zsrau * zdep  
-                    tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta        ! add the trend to the general tracer trend
-                    tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa
+                    tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + ( rnf_tmp(ji,jj) - tsn(ji,jj,jk,jp_tem) ) * zrnf
+                    tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + ( rnf_sal(ji,jj) - tsn(ji,jj,jk,jp_sal) ) * zrnf
                   END DO  
                 ENDIF  
   
-              ELSEIF( rnf(ji,jj) < 0.e0) THEN   ! for use in baltic when flow is out of domain, want no change in temp and sal  
+              ELSEIF( rnf(ji,jj) < 0.e0 ) THEN   ! for use in baltic when flow is out of domain, want no change in temp and sal  
   
                 IF( lk_vvl ) THEN  
                   ! calculate automatic adjustment to sal and temp due to dilution/concentraion effect   
-                  zata = tsn(ji,jj,1,jp_tem) * rnf(ji,jj) * zsrau * zse3t  
-                  zasa = tsn(ji,jj,1,jp_sal) * rnf(ji,jj) * zsrau * zse3t  
-                  tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + zata                  ! add the trend to the general tracer trend
-                  tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + zasa
+                  zrnf = rnf(ji,jj) * zsrau * zse3t
+                  tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + tsn(ji,jj,1,jp_tem) * zrnf 
+                  tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + tsn(ji,jj,1,jp_sal) * zrnf
                 ENDIF  
   
@@ -207 +275 @@
            END DO  
         END DO  
-
-      ELSE IF( ln_rnf ) THEN      ! Concentration dilution effect on (t,s) due to runoff without T, S and depth attributes
-
-
-        DO jj = 2, jpj
-           DO ji = fs_2, fs_jpim1   ! vector opt.
-#if ! defined key_zco
-              zse3t = 1. / fse3t(ji,jj,1)
-#endif
-              IF( lk_vvl) THEN
-                 zta =    rnf(ji,jj) * zsrau * tsn(ji,jj,1,jp_tem) * zse3t       ! & cooling/heating effect of runoff
-                 zsa =    0.e0                                                   ! No salinity concent./dilut. effect
-              ELSE
-                 zta =    0.0                                            ! temperature : heat flux 
-                 zsa =  - rnf(ji,jj) * zsrau * tsn(ji,jj,1,jp_sal) * zse3t     ! salinity :  concent./dilut. effect
-              ENDIF
-              tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + zta          ! add the trend to the general tracer trend
-              tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + zsa
-           END DO
-        END DO
- 
+        !
       ENDIF  
 

branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/istate.F90

@@ -67 +67 @@
       !! ** Purpose :   Initialization of the dynamics and tracer fields.
       !!----------------------------------------------------------------------
+      ! - ML - needed for initialization of e3t_b
+      INTEGER  ::  jk     ! dummy loop indice
 
       IF(lwp) WRITE(numout,*)
@@ -128 +130 @@
          IF( ln_zps .AND. .NOT. lk_c1d )   CALL zps_hde( nit000, jpts, tsb, gtsu, gtsv,  & ! zps: before hor. gradient
             &                                                          rhd, gru , grv  )   ! of t,s,rd at ocean bottom
-         !    
+         !   
+         ! - ML - sshn could be modified by istate_eel, so that initialization of fse3t_b is done here
+         IF( lk_vvl ) THEN
+            DO jk = 1, jpk
+               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
+            ENDDO
+         ENDIF
+         ! 
       ENDIF
       !

branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/oce.F90

@@ -36 +36 @@
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sshu_b ,  sshu_n  ,  sshu_a  !: sea surface height at u-point [m]
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sshv_b ,  sshv_n  ,  sshv_a  !: sea surface height at u-point [m]
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sshf_b ,  sshf_n  ,  sshf_a  !: sea surface height at f-point [m]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::             sshf_n             !: sea surface height at f-point [m]
    !
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   spgu, spgv                   !: horizontal surface pressure gradient

branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/opa.F90

@@ -281 +281 @@
       IF( lk_diaar5    )    CALL dia_ar5_init   ! ar5 diag
                             CALL dia_ptr_init   ! Poleward TRansports initialization
+                            CALL dia_hsb_init   ! heat content, salt content and volume budgets
                             CALL trd_mod_init   ! Mixed-layer/Vorticity/Integral constraints trends
       !

branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/step.F90

@@ -221 +221 @@
                                CALL ssh_nxt( kstp )         ! sea surface height at next time step
 
+      IF( ln_diahsb        )   CALL dia_hsb( kstp )         ! - ML - global conservation diagnostics
+
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Control and restarts

branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/step_oce.F90

@@ -85 +85 @@
    USE diahth           ! thermocline depth                (dia_hth routine)
    USE diafwb           ! freshwater budget                (dia_fwb routine)
+   USE diahsb           ! heat, salt and volume budgets    (dia_hsb routine)
    USE flo_oce          ! floats variables
    USE floats           ! floats computation               (flo_stp routine)

branches/DEV_r2106_LOCEAN2010/NEMO/TOP_SRC/PISCES/p4zmort.F90

@@ -235 +235 @@
       !!
       !! ** Method  :   Read the nampismort namelist and check the parameters
-      !!      called at the first timestep (nittrc000)
+      !!      called at the first timestep
       !!
       !! ** input   :   Namelist nampismort

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

@@ -408 +408 @@
       !!
       !! ** Method  :   Read the nampisrem namelist and check the parameters
-      !!      called at the first timestep (nittrc000)
+      !!      called at the first timestep
       !!
       !! ** input   :   Namelist nampisrem

branches/DEV_r2106_LOCEAN2010/NEMO/TOP_SRC/PISCES/p4zsink.F90

@@ -321 +321 @@
       !!
       !! ** Method  :   Read the nampiskrs namelist and check the parameters
-      !!      called at the first timestep (nittrc000)
+      !!      called at the first timestep 
       !!
       !! ** input   :   Namelist nampiskrs

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

@@ -65 +65 @@
       !!----------------------------------------------------------------------
 
-      IF( kt == nittrc000 )   CALL trc_adv_ctl          ! initialisation & control of options
+      IF( kt == nit000 )   CALL trc_adv_ctl          ! initialisation & control of options
 
-      IF( neuler == 0 .AND. kt == nittrc000 ) THEN     ! at nit000
+      IF( neuler == 0 .AND. kt == nit000 ) 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
+      ELSEIF( kt <= nit000 + nn_dttrc ) THEN          ! at nit000 or nit000+1
          r2dt(:) = 2. * rdttra(:) * FLOAT(nn_dttrc)      ! = 2 rdttra (leapfrog)
       ENDIF

branches/DEV_r2106_LOCEAN2010/NEMO/TOP_SRC/TRP/trcdmp.F90

@@ -82 +82 @@
       ! 0. Initialization (first time-step only)
       !    --------------
-      IF( kt == nittrc000 ) CALL trc_dmp_init
+      IF( kt == nit000 ) CALL trc_dmp_init
 
       IF( l_trdtrc )   ALLOCATE( ztrtrd(jpi,jpj,jpk) )   ! temporary save of trends

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

@@ -62 +62 @@
       !!----------------------------------------------------------------------
 
-      IF( kt == nittrc000 )   CALL ldf_ctl          ! initialisation & control of options
+      IF( kt == nit000 )   CALL ldf_ctl          ! initialisation & control of options
 
       IF( l_trdtrc )  THEN 

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

@@ -86 +86 @@
       !!----------------------------------------------------------------------
 
-      IF( kt == nittrc000 .AND. lwp ) THEN
+      IF( kt == nit000 .AND. lwp ) THEN
          WRITE(numout,*)
          WRITE(numout,*) 'trc_nxt : time stepping on passive tracers'
@@ -109 +109 @@
 
       ! set time step size (Euler/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)
+      IF( neuler == 0 .AND. kt ==  nit000) THEN  ;  r2dt(:) =     rdttra(:) * FLOAT( nn_dttrc )  ! at nit000             (Euler)
+      ELSEIF( kt <= nit000 + 1 )           THEN  ;  r2dt(:) = 2.* rdttra(:) * FLOAT( nn_dttrc )  ! at nit000 or nit000+1 (Leapfrog)
       ENDIF
 
@@ -120 +120 @@
 
       ! Leap-Frog + Asselin filter time stepping
-      IF( lk_vvl ) THEN   ;   CALL tra_nxt_vvl( kt, nittrc000, trb, trn, tra, jptra )      ! variable volume level (vvl) 
-      ELSE                ;   CALL tra_nxt_fix( kt, nittrc000, trb, trn, tra, jptra )      ! fixed    volume level 
+      IF( lk_vvl ) THEN   ;   CALL tra_nxt_vvl( kt, 'TRC', trb, trn, tra, jptra )      ! variable volume level (vvl) 
+      ELSE                ;   CALL tra_nxt_fix( kt,        trb, trn, tra, jptra )      ! fixed    volume level 
       ENDIF
 

branches/DEV_r2106_LOCEAN2010/NEMO/TOP_SRC/TRP/trcrad.F90

@@ -54 +54 @@
       !!----------------------------------------------------------------------
 
-      IF( kt == nittrc000 ) THEN
+      IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'trc_rad : Correct artificial negative concentrations '

branches/DEV_r2106_LOCEAN2010/NEMO/TOP_SRC/TRP/trcsbc.F90

@@ -65 +65 @@
       !! * Local declarations
       INTEGER  ::   ji, jj, jn           ! dummy loop indices
-      REAL(wp) ::   ztra, zsrau, zse3t   ! temporary scalars
+      REAL(wp) ::   zsrau, zse3t   ! temporary scalars
       REAL(wp), DIMENSION(jpi,jpj) ::   zemps  ! surface freshwater flux
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrtrd
@@ -71 +71 @@
       !!----------------------------------------------------------------------
 
-      IF( kt == nittrc000 ) THEN
+      IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'trc_sbc : Passive tracers surface boundary condition'
@@ -82 +82 @@
 #if ! defined key_offline
       ! Concentration dilution effect on tracer due to evaporation, precipitation, and river runoff
-      IF( lk_vvl ) THEN   ;   zemps(:,:) = emps(:,:) - emp(:,:)   ! volume variable
-      ELSE                ;   zemps(:,:) = emps(:,:) - rnf(:,:)   ! linear free surface 
+      IF( lk_vvl ) THEN                      ! volume variable
+         zemps(:,:) = emps(:,:) - emp(:,:)   
+!!ch         zemps(:,:) = 0.
+      ELSE                                   ! linear free surface
+         IF( ln_rnf ) THEN  ;  zemps(:,:) = emps(:,:) - rnf(:,:)   !  E-P-R
+         ELSE               ;  zemps(:,:) = emps(:,:)
+         ENDIF 
       ENDIF 
 #else
@@ -91 +96 @@
       ENDIF
 #endif
+
       ! 0. initialization
       zsrau = 1. / rau0
-#if defined key_zco
-      zse3t = 1. / e3t_0(1)
-#endif
-
       DO jn = 1, jptra
          !
          IF( l_trdtrc ) ztrtrd(:,:,:) = tra(:,:,:,jn)  ! save trends
-
+         !                                             ! add the trend to the general tracer trend
          DO jj = 2, jpj
             DO ji = fs_2, fs_jpim1   ! vector opt.
-#if ! defined key_zco 
                zse3t = 1. / fse3t(ji,jj,1)
-#endif
-               ! add the trend to the general tracer trend
-               ztra = zemps(ji,jj) *  zsrau * trn(ji,jj,1,jn) * zse3t
-               tra(ji,jj,1,jn) = tra(ji,jj,1,jn) + ztra
+               tra(ji,jj,1,jn) = tra(ji,jj,1,jn) + zemps(ji,jj) *  zsrau * trn(ji,jj,1,jn) * zse3t
             END DO
          END DO

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

@@ -61 +61 @@
       !!---------------------------------------------------------------------
 
-      IF( kt == nittrc000 )   CALL zdf_ctl          ! initialisation & control of options
+      IF( kt == nit000 )   CALL zdf_ctl          ! initialisation & control of options
 
 #if ! defined key_pisces
-      IF( neuler == 0 .AND. kt == nittrc000 ) THEN     ! at nit000
+      IF( neuler == 0 .AND. kt == nit000 ) 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
+      ELSEIF( kt <= nit000 + nn_dttrc ) THEN          ! at nit000 or nit000+1
          r2dt(:) = 2. * rdttra(:) * FLOAT(nn_dttrc)      ! = 2 rdttra (leapfrog)
       ENDIF

branches/DEV_r2106_LOCEAN2010/NEMO/TOP_SRC/TRP/trdmld_trc.F90

@@ -361 +361 @@
       !!----------------------------------------------------------------------
 
-      IF( llwarn ) THEN                                           ! warnings
-         IF(      ( nittrc000 /= nit000   ) &
-              .OR.( nn_dttrc    /= 1        )    ) THEN
-
-            WRITE(numout,*) 'Be careful, trends diags never validated'
-            STOP 'Uncomment this line to proceed'
-         ENDIF
-      ENDIF
+      IF( nn_dttrc  /= 1  ) CALL ctl_stop( " Be careful, trends diags never validated " )
 
       ! ======================================================================
@@ -422 +415 @@
       ! II.1 Set before values of vertically averages passive tracers
       ! -------------------------------------------------------------
-      IF( kt > nittrc000 ) THEN
+      IF( kt > nit000 ) THEN
          DO jn = 1, jptra
             IF( ln_trdtrc(jn) ) THEN
@@ -507 +500 @@
       tmltrd_trc(:,:,:,:) = tmltrd_trc(:,:,:,:) * rn_ucf_trc
 
-      itmod = kt - nittrc000 + 1
+      itmod = kt - nit000 + 1
       it    = kt
 
@@ -915 +908 @@
       !!----------------------------------------------------------------------
       ! ... Warnings
-      IF( llwarn ) THEN
-         IF(      ( nittrc000 /= nit000   ) &
-              .OR.( nn_dttrc    /= 1        )    ) THEN
-
-            WRITE(numout,*) 'Be careful, trends diags never validated'
-            STOP 'Uncomment this line to proceed'
-         END IF
-      END IF
+      IF( nn_dttrc  /= 1  ) CALL ctl_stop( " Be careful, trends diags never validated " )
 
       ! ======================================================================
@@ -1038 +1024 @@
 
       ! define time axis
-      itmod = kt - nittrc000 + 1
+      itmod = kt - nit000 + 1
       it    = kt
 
@@ -1318 +1304 @@
             CALL dia_nam( clhstnam, nn_trd_trc, csuff )
             CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,                                            &
-               &        1, jpi, 1, jpj, nittrc000-nn_dttrc, zjulian, rdt, nh_t(jn), nidtrd(jn), domain_id=nidom )
+               &        1, jpi, 1, jpj, nit000, zjulian, rdt, nh_t(jn), nidtrd(jn), domain_id=nidom )
       
             !-- Define the ML depth variable
@@ -1331 +1317 @@
           CALL dia_nam( clhstnam, nn_trd_trc, 'trdbio' )
           CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,                                            &
-             &             1, jpi, 1, jpj, nittrc000-nn_dttrc, zjulian, rdt, nh_tb, nidtrdbio, domain_id=nidom )
+             &             1, jpi, 1, jpj, nit000, zjulian, rdt, nh_tb, nidtrdbio, domain_id=nidom )
 #endif
 

branches/DEV_r2106_LOCEAN2010/NEMO/TOP_SRC/TRP/trdmod_trc.F90

@@ -50 +50 @@
       !!----------------------------------------------------------------------
 
-      IF( kt == nittrc000 ) THEN
+      IF( kt == nit000 ) THEN
 !         IF(lwp)WRITE(numout,*)
 !         IF(lwp)WRITE(numout,*) 'trd_mod_trc:'

branches/DEV_r2106_LOCEAN2010/NEMO/TOP_SRC/oce_trc.F90

@@ -201 +201 @@
    USE sbc_oce , ONLY :   emps       =>    emps       !: freshwater budget: concentration/dillution   [Kg/m2/s]
    USE sbc_oce , ONLY :   rnf        =>    rnf        !: river runoff   [Kg/m2/s]
+   USE sbc_oce , ONLY :   ln_rnf     =>    ln_rnf     !: runoffs / runoff mouths
    USE sbc_oce , ONLY :   fr_i       =>    fr_i       !: ice fraction (between 0 to 1)
    USE traqsr  , ONLY :   rn_abs     =>    rn_abs     !: fraction absorbed in the very near surface