Changeset 10009

Timestamp:

2018-07-29T11:23:51+02:00 (6 years ago)

Author:

gm

Message:

#1911 (ENHANCE-04): RK3 branch - step II.1 time-level dimension on ssh

Location:

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src

Files:

• ICE/iceistate.F90 (modified) (1 diff)
• ICE/iceupdate.F90 (modified) (1 diff)
• NST/agrif_oce_interp.F90 (modified) (16 diffs)
• NST/agrif_oce_update.F90 (modified) (8 diffs)
• NST/agrif_user.F90 (modified) (10 diffs)
• OCE/ASM/asmbkg.F90 (modified) (2 diffs)
• OCE/ASM/asminc.F90 (modified) (2 diffs)
• OCE/BDY/bdydta.F90 (modified) (1 diff)
• OCE/BDY/bdydyn.F90 (modified) (1 diff)
• OCE/CRS/crsfld.F90 (modified) (1 diff)
• OCE/DIA/dia25h.F90 (modified) (3 diffs)
• OCE/DIA/diaar5.F90 (modified) (5 diffs)
• OCE/DIA/diadct.F90 (modified) (7 diffs)
• OCE/DIA/diaharm.F90 (modified) (1 diff)
• OCE/DIA/diahsb.F90 (modified) (5 diffs)
• OCE/DIA/diahth.F90 (modified) (1 diff)
• OCE/DIA/diatmb.F90 (modified) (1 diff)
• OCE/DIA/diawri.F90 (modified) (3 diffs)
• OCE/DOM/domain.F90 (modified) (8 diffs)
• OCE/DOM/domvvl.F90 (modified) (22 diffs)
• OCE/DOM/iscplhsb.F90 (modified) (1 diff)
• OCE/DOM/iscplrst.F90 (modified) (11 diffs)
• OCE/DOM/istate.F90 (modified) (3 diffs)
• OCE/DOM/restart.F90 (modified) (4 diffs)
• OCE/DYN/dyncor.F90 (copied) (copied from NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynvor.F90) (14 diffs)
• OCE/DYN/dynhpg.F90 (modified) (14 diffs)
• OCE/DYN/dynnxt.F90 (modified) (3 diffs)
• OCE/DYN/dynspg.F90 (modified) (1 diff)
• OCE/DYN/dynspg_exp.F90 (modified) (2 diffs)
• OCE/DYN/dynspg_ts.F90 (modified) (10 diffs)
• OCE/DYN/sshwzv.F90 (modified) (9 diffs)
• OCE/DYN/wet_dry.F90 (modified) (5 diffs)
• OCE/IOM/iom.F90 (modified) (3 diffs)
• OCE/OBS/diaobs.F90 (modified) (2 diffs)
• OCE/OBS/obs_read_altbias.F90 (modified) (2 diffs)
• OCE/OBS/obs_readmdt.F90 (modified) (3 diffs)
• OCE/OBS/obs_sstbias.F90 (modified) (1 diff)
• OCE/SBC/sbccpl.F90 (modified) (3 diffs)
• OCE/SBC/sbcfwb.F90 (modified) (1 diff)
• OCE/SBC/sbcice_cice.F90 (modified) (2 diffs)
• OCE/SBC/sbcssm.F90 (modified) (4 diffs)
• OCE/TRA/trasbc.F90 (modified) (1 diff)
• OCE/TRD/trdpen.F90 (modified) (1 diff)
• OCE/nemogcm.F90 (modified) (1 diff)
• OCE/oce.F90 (modified) (3 diffs)
• OCE/par_oce.F90 (modified) (4 diffs)
• OCE/stpctl.F90 (modified) (2 diffs)
• OFF/dtadyn.F90 (modified) (9 diffs)
• OFF/nemogcm.F90 (modified) (1 diff)
• SAO/sao_read.F90 (modified) (5 diffs)
• TOP/oce_trc.F90 (modified) (3 diffs)
• TOP/trcsub.F90 (modified) (11 diffs)

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/iceistate.F90

@@ -415 +415 @@
       IF( ln_ice_embd ) THEN           ! embedded sea-ice: deplete the initial ssh below sea-ice area
          !
-         sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rho0
-         sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rho0
+         ssh(:,:,Nnn) = ssh(:,:,Nnn) - snwice_mass(:,:) * r1_rho0
+         ssh(:,:,Nbb) = ssh(:,:,Nbb) - snwice_mass(:,:) * r1_rho0
          !
          IF( .NOT.ln_linssh ) THEN     ! modified the now and before vertical mesh and scale factors 

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/iceupdate.F90

@@ -15 +15 @@
    !!   ice_update_tau   : update i- and j-stresses, and its modulus at the ocean surface
    !!----------------------------------------------------------------------
-   USE oce     , ONLY : sshn, sshb
    USE phycst         ! physical constants
    USE dom_oce        ! ocean domain

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/NST/agrif_oce_interp.F90

@@ -520 +520 @@
          DO jj = 1, jpj
             DO ji = 2, indx
-               ssha(ji,jj) = hbdy_w(ji-1,jj)
-            ENDDO
-         ENDDO
+               ssh(ji,jj,Naa) = hbdy_w(ji-1,jj)
+            END DO
+         END DO
       ENDIF
       !
@@ -530 +530 @@
          DO jj = 1, jpj
             DO ji = indx, nlci-1
-               ssha(ji,jj) = hbdy_e(ji-indx+1,jj)
-            ENDDO
-         ENDDO
+               ssh(ji,jj,Naa) = hbdy_e(ji-indx+1,jj)
+            END DO
+         END DO
       ENDIF
       !
@@ -540 +540 @@
          DO jj = 2, indy
             DO ji = 1, jpi
-               ssha(ji,jj) = hbdy_s(ji,jj-1)
-            ENDDO
-         ENDDO
+               ssh(ji,jj,Naa) = hbdy_s(ji,jj-1)
+            END DO
+         END DO
       ENDIF
       !
@@ -550 +550 @@
          DO jj = indy, nlcj-1
             DO ji = 1, jpi
-               ssha(ji,jj) = hbdy_n(ji,jj-indy+1)
-            ENDDO
-         ENDDO
+               ssh(ji,jj,Naa) = hbdy_n(ji,jj-indy+1)
+            END DO
+         END DO
       ENDIF
       !
@@ -576 +576 @@
             DO ji = 2, indx
                ssha_e(ji,jj) = hbdy_w(ji-1,jj)
-            ENDDO
-         ENDDO
+            END DO
+         END DO
       ENDIF
       !
@@ -586 +586 @@
             DO ji = indx, nlci-1
                ssha_e(ji,jj) = hbdy_e(ji-indx+1,jj)
-            ENDDO
-         ENDDO
+            END DO
+         END DO
       ENDIF
       !
@@ -596 +596 @@
             DO ji = 1, jpi
                ssha_e(ji,jj) = hbdy_s(ji,jj-1)
-            ENDDO
-         ENDDO
+            END DO
+         END DO
       ENDIF
       !
@@ -606 +606 @@
             DO ji = 1, jpi
                ssha_e(ji,jj) = hbdy_n(ji,jj-indy+1)
-            ENDDO
-         ENDDO
+            END DO
+         END DO
       ENDIF
       !
@@ -700 +700 @@
                   N_out = N_out + 1
                   h_out(jk) = e3t_n(iref,jref,jk)
-               ENDDO
+               END DO
                IF (N_in > 0) THEN
                   DO jn=1,jpts
                      call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out)
-                  ENDDO
+                  END DO
                ENDIF
-            ENDDO
-         ENDDO
+            END DO
+         END DO
 # else
          ptab_child(i1:i2,j1:j2,1:jpk,1:jpts) = ptab(i1:i2,j1:j2,1:jpk,1:jpts)
@@ -848 +848 @@
    END SUBROUTINE interptsn
 
+
    SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before, nb, ndir )
       !!----------------------------------------------------------------------
@@ -861 +862 @@
       !
       IF( before) THEN
-         ptab(i1:i2,j1:j2) = sshn(i1:i2,j1:j2)
+         ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Nnn)
       ELSE
          western_side  = (nb == 1).AND.(ndir == 1)
@@ -868 +869 @@
          northern_side = (nb == 2).AND.(ndir == 2)
          !! clem ghost
-         IF(western_side)  hbdy_w(1:nbghostcells,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
-         IF(eastern_side)  hbdy_e(1:nbghostcells,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
-         IF(southern_side) hbdy_s(i1:i2,1:nbghostcells) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1) 
-         IF(northern_side) hbdy_n(i1:i2,1:nbghostcells) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
+         IF(western_side )   hbdy_w(1:nbghostcells,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
+         IF(eastern_side )   hbdy_e(1:nbghostcells,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
+         IF(southern_side)   hbdy_s(i1:i2,1:nbghostcells) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1) 
+         IF(northern_side)   hbdy_n(i1:i2,1:nbghostcells) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
       ENDIF
       !
    END SUBROUTINE interpsshn
+
 
    SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir )
@@ -925 +927 @@
                   tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
                   h_in(N_in) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 
-              ENDDO
+              END DO
          
               IF (N_in == 0) THEN
@@ -937 +939 @@
                  N_out = N_out + 1
                  h_out(N_out) = e3u_a(iref,jj,jk)
-              ENDDO
+              END DO
          
               IF (N_out == 0) THEN
@@ -953 +955 @@
               ENDIF
               call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ua(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
-            ENDDO
-         ENDDO
+            END DO
+         END DO
 
 # else
@@ -1416 +1418 @@
                   N_out = N_out + 1
                   h_out(jk) = e3t_n(ji,jj,jk)
-               ENDDO
+               END DO
                IF (N_in > 0) THEN
                   CALL reconstructandremap(tabin(1:N_in),h_in,avm_k(ji,jj,1:N_out),h_out,N_in,N_out)
                ENDIF
-            ENDDO
-         ENDDO
+            END DO
+         END DO
 #else
          avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2,1)

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/NST/agrif_oce_update.F90

@@ -1015 +1015 @@
       ! 
       IF( before ) THEN
-         DO jj=j1,j2
-            DO ji=i1,i2
-               tabres(ji,jj) = sshn(ji,jj)
+         DO jj = j1, j2
+            DO ji = i1, i2
+               tabres(ji,jj) = ssh(ji,jj,Nnn)
             END DO
          END DO
@@ -1023 +1023 @@
          IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN
 !!gm         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
-            DO jj=j1,j2
-               DO ji=i1,i2
-                  sshb(ji,jj) = sshb(ji,jj) + rn_atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1)
-               END DO
-            END DO
-         ENDIF
-         !
-         DO jj=j1,j2
-            DO ji=i1,i2
-               sshn(ji,jj) = tabres(ji,jj) * tmask(ji,jj,1)
+            DO jj = j1, j2
+               DO ji = i1, i2
+                  ssh(ji,jj,Nbb) = ssh(ji,jj,Nbb) + rn_atfp * ( tabres(ji,jj) - ssh(ji,jj,Nnn) ) * tmask(ji,jj,1)
+               END DO
+            END DO
+         ENDIF
+         !
+         DO jj = j1, j2
+            DO ji = i1, i2
+               ssh(ji,jj,Nnn) = tabres(ji,jj) * tmask(ji,jj,1)
             END DO
          END DO
@@ -1038 +1038 @@
          IF ( l_1st_euler .AND. Agrif_Nb_Step() == 0 ) THEN
 !!gm         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
-            sshb(i1:i2,j1:j2)  = sshn(i1:i2,j1:j2)
+            ssh(i1:i2,j1:j2,Nbb)  = ssh(i1:i2,j1:j2,Nnn)
          ENDIF
          !
@@ -1119 +1119 @@
             DO jj=j1,j2
                zcor = rn_Dt * r1_e1e2t(i1  ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj)) 
-               sshn(i1  ,jj) = sshn(i1  ,jj) + zcor
-               IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) )   sshb(i1  ,jj) = sshb(i1  ,jj) + rn_atfp * zcor
+               ssh(i1  ,jj,Nnn) = ssh(i1  ,jj,Nnn) + zcor
+               IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) )   ssh(i1  ,jj,Nbb) = ssh(i1  ,jj,Nbb) + rn_atfp * zcor
 !!gm               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i1  ,jj) = sshb(i1  ,jj) + rn_atfp * zcor
             END DO
@@ -1127 +1127 @@
             DO jj=j1,j2
                zcor = - rn_Dt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj))
-               sshn(i2+1,jj) = sshn(i2+1,jj) + zcor
-               IF (.NOT.( lk_agrif_fstep .AND. l_1st_euler ) )   sshb(i2+1,jj) = sshb(i2+1,jj) + rn_atfp * zcor
+               ssh(i2+1,jj,Nnn) = ssh(i2+1,jj,Nnn) + zcor
+               IF (.NOT.( lk_agrif_fstep .AND. l_1st_euler ) )   ssh(i2+1,jj,Nbb) = ssh(i2+1,jj,Nbb) + rn_atfp * zcor
 !!gm               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i2+1,jj) = sshb(i2+1,jj) + rn_atfp * zcor
             END DO
@@ -1210 +1210 @@
          IF (southern_side) THEN
             DO ji=i1,i2
-               zcor = rn_Dt * r1_e1e2t(ji,j1  ) * e1v(ji,j1  ) * (vb2_b(ji,j1)-tabres(ji,j1))
-               sshn(ji,j1  ) = sshn(ji,j1  ) + zcor
-               IF ( .NOT.( lk_agrif_fstep .AND. l_euler ) )   sshb(ji,j1  ) = sshb(ji,j1) + rn_atfp * zcor
+               zcor = rn_Dt * r1_e1e2t(ji,j1  ) * e1v(ji,j1  ) * ( vb2_b(ji,j1)-tabres(ji,j1) )
+               ssh(ji,j1  ,Nnn) = ssh(ji,j1  ,Nnn) + zcor
+               IF ( .NOT.( lk_agrif_fstep .AND. l_euler ) )   ssh(ji,j1  ,Nbb) = ssh(ji,j1,Nbb) + rn_atfp * zcor
 !!gm               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j1  ) = sshb(ji,j1) + rn_atfp * zcor
             END DO
@@ -1218 +1218 @@
          IF (northern_side) THEN               
             DO ji=i1,i2
-               zcor = - rn_Dt * r1_e1e2t(ji,j2+1) * e1v(ji,j2  ) * (vb2_b(ji,j2)-tabres(ji,j2))
-               sshn(ji,j2+1) = sshn(ji,j2+1) + zcor
-               IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) )   sshb(ji,j2+1) = sshb(ji,j2+1) + rn_atfp * zcor
+               zcor = - rn_Dt * r1_e1e2t(ji,j2+1) * e1v(ji,j2  ) * ( vb2_b(ji,j2)-tabres(ji,j2) )
+               ssh(ji,j2+1,Nnn) = ssh(ji,j2+1,Nnn) + zcor
+               IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) )   ssh(ji,j2+1,Nbb) = ssh(ji,j2+1,Nbb) + rn_atfp * zcor
 !!gm               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j2+1) = sshb(ji,j2+1) + rn_atfp * zcor
             END DO
@@ -1350 +1350 @@
          ! Update e3t from ssh (z* case only)
          DO jk = 1, jpkm1
-            DO jj=j1,j2
-               DO ji=i1,i2
-                  ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + sshn(ji,jj) &
-                                     & *ssmask(ji,jj)/(ht_0(ji,jj)-1._wp + ssmask(ji,jj)))
+            DO jj = j1, j2
+               DO ji = i1, i2
+                  ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + ssh(ji,jj,Nnn) * r1_ht_0(ji,jj) *tmask(ji,jj,jk) )
                END DO
             END DO

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/NST/agrif_user.F90

@@ -190 +190 @@
    Agrif_UseSpecialValue = .TRUE.
    CALL Agrif_Bc_variable(sshn_id,calledweight=1., procname=interpsshn )
-   hbdy_w(:,:) = 0.e0 ; hbdy_e(:,:) = 0.e0 ; hbdy_n(:,:) = 0.e0 ; hbdy_s(:,:) = 0.e0
-   ssha(:,:) = 0.e0
+   hbdy_w(:,:) = 0._wp   ;   hbdy_e(:,:) = 0._wp 
+   hbdy_n(:,:) = 0._wp   ;   hbdy_s(:,:) = 0._wp
+   !
+   ssh   (:,:,Naa) = 0._wp
 
    IF ( ln_dynspg_ts ) THEN
@@ -199 +201 @@
       CALL Agrif_Bc_variable(ub2b_interp_id,calledweight=1.,procname=interpub2b)
       CALL Agrif_Bc_variable(vb2b_interp_id,calledweight=1.,procname=interpvb2b)
-      ubdy_w(:,:) = 0.e0 ; vbdy_w(:,:) = 0.e0
-      ubdy_e(:,:) = 0.e0 ; vbdy_e(:,:) = 0.e0
-      ubdy_n(:,:) = 0.e0 ; vbdy_n(:,:) = 0.e0
-      ubdy_s(:,:) = 0.e0 ; vbdy_s(:,:) = 0.e0
+      ubdy_w(:,:) = 0._wp   ;   vbdy_w(:,:) = 0._wp
+      ubdy_e(:,:) = 0._wp   ;   vbdy_e(:,:) = 0._wp
+      ubdy_n(:,:) = 0._wp   ;   vbdy_n(:,:) = 0._wp
+      ubdy_s(:,:) = 0._wp   ;   vbdy_s(:,:) = 0._wp
    ENDIF
 
    Agrif_UseSpecialValue = .FALSE. 
-   ! reset velocities to zero
-   ua(:,:,:) = 0.
-   va(:,:,:) = 0.
+   
+   ua(:,:,:) = 0._wp   ! reset velocities to zero
+   va(:,:,:) = 0._wp
 
    ! 3. Some controls
@@ -214 +216 @@
    check_namelist = .TRUE.
 
-   IF( check_namelist ) THEN 
-
-      ! Check time steps           
+   IF( check_namelist ) THEN          
+      ! Check time steps 
       IF( NINT(Agrif_Rhot()) * NINT(rn_Dt) /= Agrif_Parent(rn_Dt) ) THEN
          WRITE(cl_check1,*)  NINT(Agrif_Parent(rn_Dt))
@@ -222 +223 @@
          WRITE(cl_check3,*)  NINT(Agrif_Parent(rn_Dt)/Agrif_Rhot())
          CALL ctl_stop( 'Incompatible time step between ocean grids',   &
-               &               'parent grid value : '//cl_check1    ,   & 
-               &               'child  grid value : '//cl_check2    ,   & 
-               &               'value on child grid should be changed to : '//cl_check3 )
+            &                  'parent grid value : '//cl_check1    ,   & 
+            &                  'child  grid value : '//cl_check2    ,   & 
+            &                  'value on child grid should be changed to : '//cl_check3 )
       ENDIF
 
@@ -292 +293 @@
 END SUBROUTINE Agrif_InitValues_cont
 
+
 SUBROUTINE agrif_declare_var
    !!----------------------------------------------------------------------
@@ -450 +452 @@
 
 #if defined key_si3
+
 SUBROUTINE Agrif_InitValues_cont_ice
    !!----------------------------------------------------------------------
@@ -493 +496 @@
    !
 END SUBROUTINE Agrif_InitValues_cont_ice
+
 
 SUBROUTINE agrif_declare_var_ice
@@ -548 +552 @@
 
 END SUBROUTINE agrif_declare_var_ice
+
 #endif
 
 
 # if defined key_top
+
 SUBROUTINE Agrif_InitValues_cont_top
    !!----------------------------------------------------------------------
@@ -808 +814 @@
 
 #else
+
 SUBROUTINE Subcalledbyagrif
    !!----------------------------------------------------------------------
@@ -814 +821 @@
    WRITE(*,*) 'Impossible to be here'
 END SUBROUTINE Subcalledbyagrif
+
 #endif

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ASM/asmbkg.F90

@@ -103 +103 @@
             CALL iom_rstput( kt, nitbkg_r, inum, 'tn'     , tsn(:,:,:,jp_tem) )
             CALL iom_rstput( kt, nitbkg_r, inum, 'sn'     , tsn(:,:,:,jp_sal) )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'sshn'   , sshn              )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'sshn'   , ssh(:,:,Nnn)      )
             IF( ln_zdftke )   CALL iom_rstput( kt, nitbkg_r, inum, 'en'     , en                )
             !
@@ -138 +138 @@
             CALL iom_rstput( kt, nitdin_r, inum, 'tn'     , tsn(:,:,:,jp_tem) )
             CALL iom_rstput( kt, nitdin_r, inum, 'sn'     , tsn(:,:,:,jp_sal) )
-            CALL iom_rstput( kt, nitdin_r, inum, 'sshn'   , sshn              )
+            CALL iom_rstput( kt, nitdin_r, inum, 'sshn'   , ssh(:,:,Nnn)      )
 #if defined key_si3
             IF( nn_ice == 2 ) THEN

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ASM/asminc.F90

@@ -728 +728 @@
             IF(lwp) THEN
                WRITE(numout,*) 
-               WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', &
-                  &  kt,' with IAU weight = ', wgtiau(it)
+               WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
                WRITE(numout,*) '~~~~~~~~~~~~'
             ENDIF
@@ -755 +754 @@
          IF ( kt == nitdin_r ) THEN
             !
-            l_1st_euler = .TRUE.                         ! Force Euler forward step
-            !
-            sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:)   ! Initialize the now fields the background + increment
-            !
-            sshb(:,:) = sshn(:,:)                        ! Update before fields
-            e3t_b(:,:,:) = e3t_n(:,:,:)
+            l_1st_euler = .TRUE.                            ! Force Euler forward step
+            !
+            ssh(:,:,Nnn) = ssh_bkg(:,:) + ssh_bkginc(:,:)   ! Initialize the now fields the background + increment
+            !
+            ssh  (:,:,Nbb) = ssh  (:,:,Nnn)                 ! Update before fields
+            e3t_b(:,:,:)   = e3t_n(:,:,:)
             
 !!gm BUG :   missing the update of all other scale factors (e3u e3v e3w  etc... _n and _b) 

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/BDY/bdydta.F90

@@ -111 +111 @@
                      ii = idx_bdy(jbdy)%nbi(ib,igrd)
                      ij = idx_bdy(jbdy)%nbj(ib,igrd)
-                     dta_bdy(jbdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)         
+                     dta_bdy(jbdy)%ssh(ib) = ssh(ii,ij,Nnn) * tmask(ii,ij,1)         
                   END DO 
                ENDIF

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/BDY/bdydyn.F90

@@ -97 +97 @@
       !-------------------------------------------------------
 
-      IF( ll_dyn2d )   CALL bdy_dyn2d( kt, pua2d, pva2d, ub_b, vb_b, r1_hu_a(:,:), r1_hv_a(:,:), ssha )
+      IF( ll_dyn2d )   CALL bdy_dyn2d( kt, pua2d, pva2d, ub_b, vb_b, r1_hu_a(:,:), r1_hv_a(:,:), ssh(:,:,Naa) )
 
       IF( ll_dyn3d )   CALL bdy_dyn3d( kt )

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/CRS/crsfld.F90

@@ -219 +219 @@
       
       !  sbc fields  
-      CALL crs_dom_ope( sshn , 'VOL', 'T', tmask, sshn_crs , p_e12=e1e2t, p_e3=ze3t           , psgn=1.0 )  
-      CALL crs_dom_ope( utau , 'SUM', 'U', umask, utau_crs , p_e12=e2u  , p_surf_crs=e2u_crs  , psgn=1.0 )
-      CALL crs_dom_ope( vtau , 'SUM', 'V', vmask, vtau_crs , p_e12=e1v  , p_surf_crs=e1v_crs  , psgn=1.0 )
-      CALL crs_dom_ope( wndm , 'SUM', 'T', tmask, wndm_crs , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0 )
-      CALL crs_dom_ope( rnf  , 'MAX', 'T', tmask, rnf_crs                                     , psgn=1.0 )
-      CALL crs_dom_ope( qsr  , 'SUM', 'T', tmask, qsr_crs  , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0 )
-      CALL crs_dom_ope( emp_b, 'SUM', 'T', tmask, emp_b_crs, p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0 )
-      CALL crs_dom_ope( emp  , 'SUM', 'T', tmask, emp_crs  , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0 )
-      CALL crs_dom_ope( sfx  , 'SUM', 'T', tmask, sfx_crs  , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0 )
-      CALL crs_dom_ope( fr_i , 'SUM', 'T', tmask, fr_i_crs , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0 )
+      CALL crs_dom_ope( ssh(:,:,Nnn) , 'VOL', 'T', tmask, sshn_crs , p_e12=e1e2t, p_e3=ze3t           , psgn=1.0 ) 
+      CALL crs_dom_ope( utau         , 'SUM', 'U', umask, utau_crs , p_e12=e2u  , p_surf_crs=e2u_crs  , psgn=1.0 )
+      CALL crs_dom_ope( vtau         , 'SUM', 'V', vmask, vtau_crs , p_e12=e1v  , p_surf_crs=e1v_crs  , psgn=1.0 )
+      CALL crs_dom_ope( wndm         , 'SUM', 'T', tmask, wndm_crs , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0 )
+      CALL crs_dom_ope( rnf          , 'MAX', 'T', tmask, rnf_crs                                     , psgn=1.0 )
+      CALL crs_dom_ope( qsr          , 'SUM', 'T', tmask, qsr_crs  , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0 )
+      CALL crs_dom_ope( emp_b        , 'SUM', 'T', tmask, emp_b_crs, p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0 )
+      CALL crs_dom_ope( emp          , 'SUM', 'T', tmask, emp_crs  , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0 )
+      CALL crs_dom_ope( sfx          , 'SUM', 'T', tmask, sfx_crs  , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0 )
+      CALL crs_dom_ope( fr_i         , 'SUM', 'T', tmask, fr_i_crs , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0 )
 
       CALL iom_put( "ssh"      , sshn_crs )   ! ssh output 

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/dia25h.F90

@@ -94 +94 @@
       ! ------------------------- !
       cnt_25h = 1  ! sets the first value of sum at timestep 1 (note - should strictly be at timestep zero so before values used where possible) 
-      tn_25h  (:,:,:) = tsb (:,:,:,jp_tem)
-      sn_25h  (:,:,:) = tsb (:,:,:,jp_sal)
-      sshn_25h(:,:)   = sshb(:,:)
-      un_25h  (:,:,:) = ub  (:,:,:)
-      vn_25h  (:,:,:) = vb  (:,:,:)
-      wn_25h  (:,:,:) = wn  (:,:,:)
-      avt_25h (:,:,:) = avt (:,:,:)
-      avm_25h (:,:,:) = avm (:,:,:)
+      tn_25h  (:,:,:) = tsb(:,:,:,jp_tem)
+      sn_25h  (:,:,:) = tsb(:,:,:,jp_sal)
+      sshn_25h(:,:)   = ssh(:,:,Nbb)
+      un_25h  (:,:,:) = ub (:,:,:)
+      vn_25h  (:,:,:) = vb (:,:,:)
+      wn_25h  (:,:,:) = wn (:,:,:)
+      avt_25h (:,:,:) = avt(:,:,:)
+      avm_25h (:,:,:) = avm(:,:,:)
       IF( ln_zdftke ) THEN
          en_25h(:,:,:) = en(:,:,:)
@@ -156 +156 @@
          ENDIF
 
-         tn_25h  (:,:,:)     = tn_25h  (:,:,:) + tsn (:,:,:,jp_tem)
-         sn_25h  (:,:,:)     = sn_25h  (:,:,:) + tsn (:,:,:,jp_sal)
-         sshn_25h(:,:)       = sshn_25h(:,:)   + sshn(:,:)
-         un_25h  (:,:,:)     = un_25h  (:,:,:) + un  (:,:,:)
-         vn_25h  (:,:,:)     = vn_25h  (:,:,:) + vn  (:,:,:)
-         wn_25h  (:,:,:)     = wn_25h  (:,:,:) + wn  (:,:,:)
-         avt_25h (:,:,:)     = avt_25h (:,:,:) + avt (:,:,:)
-         avm_25h (:,:,:)     = avm_25h (:,:,:) + avm (:,:,:)
+         tn_25h  (:,:,:)     = tn_25h  (:,:,:) + tsn(:,:,:,jp_tem)
+         sn_25h  (:,:,:)     = sn_25h  (:,:,:) + tsn(:,:,:,jp_sal)
+         sshn_25h(:,:)       = sshn_25h(:,:)   + ssh(:,:,Nnn)
+         un_25h  (:,:,:)     = un_25h  (:,:,:) + un (:,:,:)
+         vn_25h  (:,:,:)     = vn_25h  (:,:,:) + vn (:,:,:)
+         wn_25h  (:,:,:)     = wn_25h  (:,:,:) + wn (:,:,:)
+         avt_25h (:,:,:)     = avt_25h (:,:,:) + avt(:,:,:)
+         avm_25h (:,:,:)     = avm_25h (:,:,:) + avm(:,:,:)
          IF( ln_zdftke ) THEN
             en_25h(:,:,:)    = en_25h  (:,:,:) + en(:,:,:)
@@ -206 +206 @@
          zmdi=1.e+20 !missing data indicator for masking
          ! write tracers (instantaneous)
-         zw3d(:,:,:) = tn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
-         CALL iom_put("temper25h", zw3d)   ! potential temperature
-         zw3d(:,:,:) = sn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
-         CALL iom_put( "salin25h", zw3d  )   ! salinity
-         zw2d(:,:) = sshn_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1))
-         CALL iom_put( "ssh25h", zw2d )   ! sea surface 
+         zw3d(:,:,:) = tn_25h  (:,:,:)* tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))   ;   CALL iom_put("temper25h"  , zw3d )   ! potential temperature
+         zw3d(:,:,:) = sn_25h  (:,:,:)* tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))   ;   CALL iom_put( "salin25h"  , zw3d )   ! salinity
+         zw2d(:,:)   = sshn_25h(:,:)  *ssmask(:,:)   + zmdi*(1.0-tmask(:,:,1))   ;   CALL iom_put( "ssh25h"    , zw2d )   ! sea surface 
          ! Write velocities (instantaneous)
-         zw3d(:,:,:) = un_25h(:,:,:)*umask(:,:,:) + zmdi*(1.0-umask(:,:,:))
-         CALL iom_put("vozocrtx25h", zw3d)    ! i-current
-         zw3d(:,:,:) = vn_25h(:,:,:)*vmask(:,:,:) + zmdi*(1.0-vmask(:,:,:))
-         CALL iom_put("vomecrty25h", zw3d  )   ! j-current
-         zw3d(:,:,:) = wn_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
-         CALL iom_put("vomecrtz25h", zw3d )   ! k-current
+         zw3d(:,:,:) = un_25h  (:,:,:)* umask(:,:,:) + zmdi*(1.0-umask(:,:,:))   ;   CALL iom_put("vozocrtx25h", zw3d )   ! i-current
+         zw3d(:,:,:) = vn_25h  (:,:,:)* vmask(:,:,:) + zmdi*(1.0-vmask(:,:,:))   ;   CALL iom_put("vomecrty25h", zw3d )   ! j-current
+         zw3d(:,:,:) = wn_25h  (:,:,:)* wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))   ;   CALL iom_put("vomecrtz25h", zw3d )   ! k-current
          ! Write vertical physics
-         zw3d(:,:,:) = avt_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
-         CALL iom_put("avt25h", zw3d )   ! diffusivity
-         zw3d(:,:,:) = avm_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
-         CALL iom_put("avm25h", zw3d)   ! viscosity
-         IF( ln_zdftke ) THEN
-            zw3d(:,:,:) = en_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
-            CALL iom_put("tke25h", zw3d)   ! tke
-         ENDIF
-         IF( ln_zdfgls ) THEN
-            zw3d(:,:,:) = en_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
-            CALL iom_put("tke25h", zw3d)   ! tke
-            zw3d(:,:,:) = rmxln_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
-            CALL iom_put( "mxln25h",zw3d)
+         zw3d(:,:,:) = avt_25h (:,:,:)* wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))   ;   CALL iom_put("avt25h"     , zw3d )   ! diffusivity
+         zw3d(:,:,:) = avm_25h (:,:,:)* wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))   ;   CALL iom_put("avm25h"     , zw3d )   ! viscosity
+         IF( ln_zdftke ) THEN
+            zw3d(:,:,:) = en_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))   ;   CALL iom_put("tke25h"     , zw3d )   ! tke
+         ENDIF
+         IF( ln_zdfgls ) THEN
+            zw3d(:,:,:) = en_25h   (:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))   ;   CALL iom_put("tke25h"  , zw3d )   ! tke
+            zw3d(:,:,:) = rmxln_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))   ;   CALL iom_put( "mxln25h", zw3d )
          ENDIF
          !
          ! After the write reset the values to cnt=1 and sum values equal current value 
-         tn_25h  (:,:,:) = tsn (:,:,:,jp_tem)
-         sn_25h  (:,:,:) = tsn (:,:,:,jp_sal)
-         sshn_25h(:,:)   = sshn(:,:)
-         un_25h  (:,:,:) = un  (:,:,:)
-         vn_25h  (:,:,:) = vn  (:,:,:)
-         wn_25h  (:,:,:) = wn  (:,:,:)
-         avt_25h (:,:,:) = avt (:,:,:)
-         avm_25h (:,:,:) = avm (:,:,:)
+         tn_25h  (:,:,:) = tsn(:,:,:,jp_tem)
+         sn_25h  (:,:,:) = tsn(:,:,:,jp_sal)
+         sshn_25h(:,:)   = ssh(:,:,Nnn)
+         un_25h  (:,:,:) = un (:,:,:)
+         vn_25h  (:,:,:) = vn (:,:,:)
+         wn_25h  (:,:,:) = wn (:,:,:)
+         avt_25h (:,:,:) = avt(:,:,:)
+         avm_25h (:,:,:) = avm(:,:,:)
          IF( ln_zdftke ) THEN
             en_25h(:,:,:) = en(:,:,:)

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/diaar5.F90

@@ -89 +89 @@
          ALLOCATE( zrhd(jpi,jpj,jpk) , zrhop(jpi,jpj,jpk) )
          ALLOCATE( ztsn(jpi,jpj,jpk,jpts) )
-         zarea_ssh(:,:) = area(:,:) * sshn(:,:)
+         zarea_ssh(:,:) = area(:,:) * ssh(:,:,Nnn)
       ENDIF
       !
@@ -100 +100 @@
          CALL iom_put( 'voltot', zvol               )
          CALL iom_put( 'sshtot', zvolssh / area_tot )
-         CALL iom_put( 'sshdyn', sshn(:,:) - (zvolssh / area_tot) )
+         CALL iom_put( 'sshdyn', ssh(:,:,Nnn) - (zvolssh / area_tot) )
          !
       ENDIF
@@ -118 +118 @@
                DO ji = 1, jpi
                   DO jj = 1, jpj
-                     zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
+                     zbotpres(ji,jj) = zbotpres(ji,jj) + ssh(ji,jj,Nnn) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
                   END DO
                END DO
             ELSE
-               zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
+               zbotpres(:,:) = zbotpres(:,:) + ssh(:,:,Nnn) * zrhd(:,:,1)
             END IF
 !!gm
@@ -147 +147 @@
                DO ji = 1,jpi
                   DO jj = 1,jpj
-                     zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
+                     zbotpres(ji,jj) = zbotpres(ji,jj) + ssh(ji,jj,Nnn) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
                   END DO
                END DO
             ELSE
-               zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
+               zbotpres(:,:) = zbotpres(:,:) + ssh(:,:,Nnn) * zrhd(:,:,1)
             END IF
          END IF
@@ -162 +162 @@
          !                                         ! ocean bottom pressure
          zztmp = rho0 * grav * 1.e-4_wp               ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa
-         zbotpres(:,:) = zztmp * ( zbotpres(:,:) + sshn(:,:) + thick0(:,:) )
+         zbotpres(:,:) = zztmp * ( zbotpres(:,:) + ssh(:,:,Nnn) + thick0(:,:) )
          CALL iom_put( 'botpres', zbotpres )
          !

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/diadct.F90

@@ -491 +491 @@
      !!
      !----------------------------------------------------------------------------
-     !! * arguments
      TYPE(SECTION),INTENT(INOUT) :: sec
      CHARACTER(len=1),INTENT(IN) :: cdind   ! = 'I'/'J'
      CHARACTER(len=8),INTENT(IN) :: cdextr  ! = 'top_list'/'bot_list'
      LOGICAL,INTENT(IN)          :: ld_debug                     
-
-     !! * Local variables
+     !
      INTEGER :: iextr         ,& !extremity of listpoint that we verify
                 iind          ,& !coord     of listpoint that we verify
@@ -588 +586 @@
      REAL(wp)::   zumid, zvmid, zumid_ice, zvmid_ice   ! U/V ocean & ice velocity on a cell segment 
      REAL(wp)::   zTnorm                               ! transport of velocity through one cell's sides 
-     REAL(wp)::   ztn, zsn, zrhoi, zrhop, zsshn, zdep  ! temperature/salinity/potential density/ssh/depth at u/v point
+     REAL(wp)::   ztn, zsn, zrhoi, zrhop, zssh, zdep   ! temperature/salinity/potential density/ssh/depth at u/v point
      TYPE(POINT_SECTION) ::   k
       !!--------------------------------------------------------
@@ -680 +678 @@
                   zrhop = interp(k%I,k%J,jk,'V',rhop) 
                   zrhoi = interp(k%I,k%J,jk,'V',rhd*rho0+rho0) 
-                  zsshn =  0.5*( sshn(k%I,k%J) + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1) 
+                  zssh  =  0.5*( ssh(k%I,k%J,Nnn) + ssh(k%I,k%J+1,Nnn)    ) * vmask(k%I,k%J,1) 
                CASE(2,3) 
                   ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 
@@ -686 +684 @@
                   zrhop = interp(k%I,k%J,jk,'U',rhop) 
                   zrhoi = interp(k%I,k%J,jk,'U',rhd*rho0+rho0) 
-                  zsshn =  0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)  
+                  zssh  =  0.5*( ssh(k%I,k%J,Nnn) + ssh(k%I+1,k%J,Nnn)    ) * umask(k%I,k%J,1)  
                END SELECT 
                !
@@ -706 +704 @@
 
 !!gm  THIS is WRONG  no transport due to ssh in linear free surface case !!!!!
-               IF( ln_linssh ) THEN              !add transport due to free surface 
-                  IF( jk==1 ) THEN 
-                     zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk)   & 
-                        &            + zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk) 
-                  ENDIF 
-               ENDIF
+!               IF( ln_linssh ) THEN              !add transport due to free surface 
+!                  IF( jk==1 ) THEN 
+!                     zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zssh * umask(k%I,k%J,jk)   & 
+!                        &            + zvmid* e1v(k%I,k%J) * zssh * vmask(k%I,k%J,jk) 
+!                  ENDIF 
+!               ENDIF
 !!gm end
               !COMPUTE TRANSPORT  
@@ -792 +790 @@
      TYPE(POINT_SECTION) :: k 
      INTEGER  :: jk,jseg,jclass                        ! dummy variables for looping on level/segment/classes  
-     REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zdep ! temperature/salinity/ssh/potential density /depth at u/v point 
+     REAL(wp) :: ztn, zsn, zrhoi, zrhop, zssh, zdep ! temperature/salinity/ssh/potential density /depth at u/v point 
      !!------------------------------------------------------------- 
  
@@ -858 +856 @@
                  zrhop = interp(k%I,k%J,jk,'U',rhop) 
                  zrhoi = interp(k%I,k%J,jk,'U',rhd*rho0+rho0) 
-                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)  
+                 zssh  =  0.5*(  ssh(k%I,k%J,Nnn) + ssh(k%I+1,k%J,Nnn)  ) * umask(k%I,k%J,1)  
               END SELECT 
  

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/diaharm.F90

@@ -193 +193 @@
                   DO ji = 1,jpi
                      ! Elevation
-                     ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*ssmask (ji,jj)        
-                     ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*ssumask(ji,jj)
-                     ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj)
+                     ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*ssh (ji,jj,Nnn)*ssmask (ji,jj)        
+                     ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)    *ssumask(ji,jj)
+                     ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)    *ssvmask(ji,jj)
                   END DO
                END DO

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/diahsb.F90

@@ -135 +135 @@
 
       !                    ! volume variation (calculated with ssh)
-      zdiff_v1 = glob_sum_full( surf(:,:)*sshn(:,:) - surf_ini(:,:)*ssh_ini(:,:) )
+      zdiff_v1 = glob_sum_full( surf(:,:)*ssh(:,:,Nnn) - surf_ini(:,:)*ssh_ini(:,:) )
 
       !                    ! heat & salt content variation (associated with ssh)
@@ -142 +142 @@
             DO ji = 1, jpi
                DO jj = 1, jpj
-                  z2d0(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj) - ssh_hc_loc_ini(ji,jj) ) 
-                  z2d1(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj) - ssh_sc_loc_ini(ji,jj) ) 
+                  z2d0(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_tem) * ssh(ji,jj,Nnn) - ssh_hc_loc_ini(ji,jj) )
+                  z2d1(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_sal) * ssh(ji,jj,Nnn) - ssh_sc_loc_ini(ji,jj) )
                END DO
             END DO
          ELSE                          ! no under ice-shelf seas
-            z2d0(:,:) = surf(:,:) * ( tsn(:,:,1,jp_tem) * sshn(:,:) - ssh_hc_loc_ini(:,:) ) 
-            z2d1(:,:) = surf(:,:) * ( tsn(:,:,1,jp_sal) * sshn(:,:) - ssh_sc_loc_ini(:,:) ) 
+            z2d0(:,:) = surf(:,:) * ( tsn(:,:,1,jp_tem) * ssh(:,:,Nnn) - ssh_hc_loc_ini(:,:) ) 
+            z2d1(:,:) = surf(:,:) * ( tsn(:,:,1,jp_sal) * ssh(:,:,Nnn) - ssh_sc_loc_ini(:,:) ) 
          END IF
          z_ssh_hc = glob_sum_full( z2d0 ) 
@@ -191 +191 @@
 !!gm to be added ?
 !      IF( ln_linssh ) THEN            ! fixed volume, add the ssh contribution
-!        zvol_tot = zvol_tot + glob_sum( surf(:,:) * sshn(:,:) )
+!        zvol_tot = zvol_tot + glob_sum( surf(:,:) * ssh(:,:,Nnn) )
 !      ENDIF
 !!gm end
@@ -281 +281 @@
             IF(lwp) WRITE(numout,*)
             surf_ini(:,:) = e1e2t(:,:) * tmask_i(:,:)         ! initial ocean surface
-            ssh_ini(:,:) = sshn(:,:)                          ! initial ssh
+            ssh_ini (:,:) = ssh(:,:,Nnn)                      ! initial ssh
             DO jk = 1, jpk
               ! if ice sheet/oceqn coupling, need to mask ini variables here (mask could change at the next NEMO instance).
@@ -295 +295 @@
                   DO ji = 1, jpi
                      DO jj = 1, jpj
-                        ssh_hc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj)   ! initial heat content in ssh
-                        ssh_sc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj)   ! initial salt content in ssh
+                        ssh_hc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) * ssh(ji,jj,Nnn)   ! initial heat content in ssh
+                        ssh_sc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) * ssh(ji,jj,Nnn)   ! initial salt content in ssh
                      END DO
                    END DO
                 ELSE
-                  ssh_hc_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:)   ! initial heat content in ssh
-                  ssh_sc_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:)   ! initial salt content in ssh
+                  ssh_hc_loc_ini(:,:) = tsn(:,:,1,jp_tem) * ssh(:,:,Nnn)   ! initial heat content in ssh
+                  ssh_sc_loc_ini(:,:) = tsn(:,:,1,jp_sal) * ssh(:,:,Nnn)   ! initial salt content in ssh
                END IF
                frc_wn_t = 0._wp                                       ! initial heat content misfit due to free surface

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/diahth.F90

@@ -311 +311 @@
       END DO
       ! surface boundary condition
-      IF( ln_linssh ) THEN   ;   zthick(:,:) = sshn(:,:)   ;   htc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)  
-      ELSE                   ;   zthick(:,:) = 0._wp       ;   htc3(:,:) = 0._wp                                   
+      IF( ln_linssh ) THEN   ;   zthick(:,:) = ssh(:,:,Nnn)   ;   htc3(:,:) = tsn(:,:,1,jp_tem) * ssh(:,:,Nnn) * tmask(:,:,1)
+      ELSE                   ;   zthick(:,:) = 0._wp          ;   htc3(:,:) = 0._wp                                   
       ENDIF
       ! integration down to ilevel

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/diatmb.F90

@@ -108 +108 @@
       CALL dia_calctmb( tsn(:,:,:,jp_tem), zwtmb )
       !ssh already output but here we output it masked
-      CALL iom_put( "sshnmasked", sshn(:,:)*tmask(:,:,1) + zmdi*(1.0 - tmask(:,:,1)) )
+      CALL iom_put( "sshnmasked", ssh(:,:,Nnn)*ssmask(:,:) + zmdi*(1._wp - ssmask(:,:)) )
       CALL iom_put( "top_temp"  , zwtmb(:,:,1) )    ! tmb Temperature
       CALL iom_put( "mid_temp"  , zwtmb(:,:,2) )    ! tmb Temperature

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/diawri.F90

@@ -136 +136 @@
 
       IF( ll_wd ) THEN
-         CALL iom_put( "ssh" , (sshn+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying)
+         CALL iom_put( "ssh" , ( ssh(:,:,Nnn)+ssh_ref)*ssmask(:,:) )   ! sea surface height (brought back to the reference used for wetting and drying)
       ELSE
-         CALL iom_put( "ssh" , sshn )              ! sea surface height
+         CALL iom_put( "ssh" ,   ssh(:,:,Nnn) )              ! sea surface height
       ENDIF
 
       IF( iom_use("wetdep") )   &                  ! wet depth
-         CALL iom_put( "wetdep" , ht_0(:,:) + sshn(:,:) )
+         CALL iom_put( "wetdep" , ht_0(:,:) + ssh(:,:,Nnn) )
       
       CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature
@@ -771 +771 @@
          CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
       ENDIF
-      CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height
+      CALL histwrite( nid_T, "sossheig", it, ssh(:,:,Nnn)  , ndim_hT, ndex_hT )   ! sea surface height
       CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
       CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs
@@ -993 +993 @@
       CALL histwrite( id_i, "votemper", kt, tsn(:,:,:,jp_tem), jpi*jpj*jpk, idex )    ! now temperature
       CALL histwrite( id_i, "vosaline", kt, tsn(:,:,:,jp_sal), jpi*jpj*jpk, idex )    ! now salinity
-      CALL histwrite( id_i, "sossheig", kt, sshn             , jpi*jpj    , idex )    ! sea surface height
+      CALL histwrite( id_i, "sossheig", kt, ssh(:,:,Nnn)     , jpi*jpj    , idex )    ! sea surface height
       CALL histwrite( id_i, "vozocrtx", kt, un               , jpi*jpj*jpk, idex )    ! now i-velocity
       CALL histwrite( id_i, "vomecrty", kt, vn               , jpi*jpj*jpk, idex )    ! now j-velocity

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domain.F90

@@ -15 +15 @@
    !!            3.7  !  2015-11  (G. Madec, A. Coward)  time varying zgr by default
    !!            4.0  !  2016-10  (G. Madec, S. Flavoni)  domain configuration / user defined interface
+   !!            5.0  !  2018-07  (G. Madec)  RK3: add a time-level dimension to time-varying variables
    !!----------------------------------------------------------------------
    
@@ -48 +49 @@
    PRIVATE
 
+   PUBLIC   dom_nam      ! called by nemogcm.F90
    PUBLIC   dom_init     ! called by nemogcm.F90
    PUBLIC   domain_cfg   ! called by nemogcm.F90
@@ -118 +120 @@
       !
       CALL dom_glo                     ! global domain versus local domain
-      CALL dom_nam                     ! read namelist ( namrun, namdom )
-      !
-      IF( lwxios ) THEN
-!define names for restart write and set core output (restart.F90)
+      !
+      IF( lwxios ) THEN                ! define names for restart write and set core output (restart.F90)
          CALL iom_set_rst_vars(rst_wfields)
          CALL iom_set_rstw_core(cdstr)
       ENDIF
-!reset namelist for SAS
-      IF(cdstr == 'SAS') THEN
+      IF( cdstr == 'SAS' ) THEN        ! reset namelist for SAS
          IF(lrxios) THEN
                IF(lwp) write(numout,*) 'Disable reading restart file using XIOS for SAS'
@@ -184 +183 @@
          r1_hu_b = r1_hu_0   ;   r1_hu_n = r1_hu_0   ; r1_hu_a = r1_hu_0   ! 1 / water column
          r1_hv_b = r1_hv_0   ;   r1_hv_n = r1_hv_0   ; r1_hv_a = r1_hv_0   !       thickness
-         !
          !
       ELSE                       != time varying : initialize before/now/after variables
@@ -307 +305 @@
       ENDIF
       !
-      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
+      REWIND( numnam_ref )                ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
       READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
 903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdom in reference namelist', lwp )
-      REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
+      REWIND( numnam_cfg )                ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
       READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
 904   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
@@ -316 +314 @@
       !
       IF(lwp) THEN
-         WRITE(numout,*)
          WRITE(numout,*) '   Namelist : namdom   ---   space & time domain'
          WRITE(numout,*) '      3rd order Runge-Kutta scheme            ln_RK3      = ', ln_RK3
@@ -337 +334 @@
       ENDIF
       !
-      REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run
+      !                                   ! Set the number of time-level
+      !
+      IF    ( ln_RK3 ) THEN                  ! two-level time-stepping scheme
+         Nt = 2                                    ! number of time-level
+         Nbb = 1     ! before=now                  ! time indexes
+         Nnn = 1     ! now   (n  )
+         Naa = 2     ! after (n+1)
+         !
+      ELSEIF( ln_MLF ) THEN                  ! three-level time-stepping scheme
+         Nt = 3                                    ! number of time-level
+         Nbb = 1     ! before (n-1)
+         Nnn = 2     ! now    (n  )
+         Naa = 3     ! after  (n+1)
+      ENDIF
+      !
+      REWIND( numnam_ref )                ! Namelist namrun in reference namelist : Parameters of the run
       READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namrun in reference namelist', lwp )
-      REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run
+      REWIND( numnam_cfg )                ! Namelist namrun in configuration namelist : Parameters of the run
       READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp )
       IF(lwm) WRITE ( numond, namrun )
       !
-      IF(lwp) THEN                  ! control print
+      IF(lwp) THEN                        ! control print
          WRITE(numout,*)
          WRITE(numout,*) '   Namelist : namrun   ---   run parameters'
@@ -383 +395 @@
       ENDIF
 
-      cexper = cn_exp         ! conversion DOCTOR names into model names (this should disappear soon)
+      cexper = cn_exp         ! conversion from namelist names to model names (this should disappear)
       nrstdt = nn_rstctl
       nit000 = nn_it000

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domvvl.F90

@@ -135 +135 @@
       IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
       !
-      !                    ! Read or initialize e3t_(b/n), ssh(b/n)
+      !                    ! Read or initialize ssh(Nbb) & ssh(Nnn)
       CALL dom_vvl_rst( nit000, 'READ' )
       !
@@ -142 +142 @@
       !                          !* BEFORE fields : 
       CALL ssh2e3_before               ! set:      hu , hv , r1_hu, r1_hv 
-      !                                    !  e3t, e3w, e3u, e3uw, e3v, e3vw
-      !
-      !                                ! set one for all last level to the e3._0 value
+      !                                !      e3t, e3w, e3u, e3uw, e3v, e3vw        (from 1 to jpkm1)
+      !
+      !                                ! set jpk level one to the e3._0 values
       e3t_b(:,:,jpk) = e3t_0(:,:,jpk)  ;   e3u_b(:,:,jpk) =  e3w_0(:,:,jpk)  ;   e3v_b(:,:,jpk) =  e3v_0(:,:,jpk)
       e3w_b(:,:,jpk) = e3w_0(:,:,jpk)  ;  e3uw_b(:,:,jpk) = e3uw_0(:,:,jpk)  ;  e3vw_b(:,:,jpk) = e3vw_0(:,:,jpk)
@@ -150 +150 @@
       !                          !* NOW fields : 
       CALL ssh2e3_now                  ! set: ht , hu , hv , r1_hu, r1_hv
-      !                                !      e3t, e3w, e3u, e3uw, e3v, e3vw, e3f
+      !                                !      e3t, e3w, e3u, e3uw, e3v, e3vw, e3f   (from 1 to jpkm1)
       !                                !      gdept_n, gdepw_n, gde3w_n
+!!gm issue?   gdept_n, gdepw_n, gde3w_n never defined at jpk
       !
       !                                ! set one for all last level to the e3._0 value
@@ -230 +231 @@
       !                                   !==  after ssh  ==!  (u- and v-points)
       DO jj = 2, jpjm1   ;   DO ji = 2, jpim1
-         zsshu_h(ji,jj) = 0.5_wp * ( ssha(ji,jj) + ssha(ji+1,jj) ) * ssumask(ji,jj)
-         zsshv_h(ji,jj) = 0.5_wp * ( ssha(ji,jj) + ssha(ji,jj+1) ) * ssvmask(ji,jj)
+         zsshu_h(ji,jj) = 0.5_wp * ( ssh(ji,jj,Naa) + ssh(ji+1,jj,Naa) ) * ssumask(ji,jj)
+         zsshv_h(ji,jj) = 0.5_wp * ( ssh(ji,jj,Naa) + ssh(ji,jj+1,Naa) ) * ssvmask(ji,jj)
       END DO             ;   END DO      
       CALL lbc_lnk_multi( zsshu_h(:,:), 'U', 1._wp , zsshv_h(:,:), 'V', 1._wp )
@@ -242 +243 @@
       !
       !                                   !==  after scale factors  ==!  (e3t , e3u , e3v)
-      zssht_h(:,:) = ssha   (:,:) * r1_ht_0(:,:)           ! t-point
-      zsshu_h(:,:) = zsshu_h(:,:) * r1_hu_0(:,:)           ! u-point
-      zsshv_h(:,:) = zsshv_h(:,:) * r1_hv_0(:,:)           ! v-point
+      zssht_h(:,:) = ssh    (:,:,Naa) * r1_ht_0(:,:)     ! t-point
+      zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)     ! u-point
+      zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)     ! v-point
       DO jk = 1, jpkm1
          e3t_a(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
@@ -318 +319 @@
          e3v_n  (:,:,jk) = e3v_a  (:,:,jk)
       END DO
-      ht_n(:,:) = ht_0(:,:) + sshn(:,:)            ! ocean thickness
+      ht_n(:,:) = ht_0(:,:) + ssh(:,:,Nnn)            ! ocean thickness
       !
       hu_n(:,:) = hu_a(:,:)   ;   r1_hu_n(:,:) = r1_hu_a(:,:)
@@ -326 +327 @@
       !                                            !* ssh at u- and v-points)
       DO jj = 2, jpjm1   ;   DO ji = 2, jpim1
-         zsshu_h(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji+1,jj  ) ) * ssumask(ji,jj)
-         zsshv_h(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji  ,jj+1) ) * ssvmask(ji,jj)
+         zsshu_h(ji,jj) = 0.5_wp  * ( ssh(ji,jj,Nbb) + ssh(ji+1,jj  ,Nbb) ) * ssumask(ji,jj)
+         zsshv_h(ji,jj) = 0.5_wp  * ( ssh(ji,jj,Nbb) + ssh(ji  ,jj+1,Nbb) ) * ssvmask(ji,jj)
       END DO             ;   END DO      
       CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp )
       !
       !                                            !*  e3w_b , e3uw_b , e3vw_b
-      zssht_h(:,:) = sshb (:,:) * r1_ht_0(:,:)           ! w-point
-      zsshu_h(:,:) = zsshu_h(:,:) * r1_hu_0(:,:)           ! uw-point
-      zsshv_h(:,:) = zsshv_h(:,:) * r1_hv_0(:,:)           ! vw-point
+      zssht_h(:,:) = ssh    (:,:,Nbb) * r1_ht_0(:,:)     ! w-point
+      zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)     ! uw-point
+      zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)     ! vw-point
       DO jk = 1, jpkm1
-          e3w_b(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
+          e3w_b(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk) , tmask(:,:,jk+1) ) )
          e3uw_b(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) )
          e3vw_b(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) )
@@ -344 +345 @@
       !                                            !* ssh at u- and v-points)
       DO jj = 1, jpjm1   ;   DO ji = 1, jpim1            ! start from 1 for f-point
-         zsshu_h(ji,jj) = 0.50_wp * ( sshn(ji  ,jj) + sshn(ji+1,jj  ) ) * ssumask(ji,jj)
-         zsshv_h(ji,jj) = 0.50_wp * ( sshn(ji  ,jj) + sshn(ji  ,jj+1) ) * ssvmask(ji,jj)
-         zsshf_h(ji,jj) = 0.25_wp * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)   & 
-            &                       + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) * ssfmask(ji,jj)
+         zsshu_h(ji,jj) = 0.50_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji+1,jj  ,Nnn) ) * ssumask(ji,jj)
+         zsshv_h(ji,jj) = 0.50_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn) ) * ssvmask(ji,jj)
+         zsshf_h(ji,jj) = 0.25_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn)   & 
+            &                       + ssh(ji+1,jj,Nnn) + ssh(ji+1,jj+1,Nnn) ) * ssfmask(ji,jj)
       END DO             ;   END DO      
       CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp , zsshf_h(:,:),'F', 1._wp )      
       !
       !                                            !* e3w_n , e3uw_n , e3vw_n, e3f_n 
-      zssht_h(:,:) = sshn   (:,:) * r1_ht_0(:,:)           ! t- & w-point
-      zsshu_h(:,:) = zsshu_h(:,:) * r1_hu_0(:,:)           ! uw-point
-      zsshv_h(:,:) = zsshv_h(:,:) * r1_hv_0(:,:)           ! vw-point
-      zsshf_h(:,:) = zsshf_h(:,:) * r1_hf_0(:,:)           ! f-point
+      zssht_h(:,:) = ssh    (:,:,Nnn) * r1_ht_0(:,:)     ! t- & w-point
+      zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)     ! uw-point
+      zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)     ! vw-point
+      zsshf_h(:,:) = zsshf_h(:,:)     * r1_hf_0(:,:)     ! f-point
       DO jk = 1, jpkm1
-          e3w_n(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) *  tmask(:,:,jk) )
+          e3w_n(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk) , tmask(:,:,jk+1) ) )
          e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * wumask(:,:,jk) )
          e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * wvmask(:,:,jk) )
@@ -363 +364 @@
       END DO      
       ! 
-      zssht_h(:,:) = 1._wp + sshn (:,:) * r1_ht_0(:,:)   ! t-point
+      zssht_h(:,:) = 1._wp + zssht_h(:,:)               ! t-point
       !
       IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness 
@@ -369 +370 @@
             gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:)
             gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:)
-            gde3w_n(:,:,jk) =   gdept_n(:,:,jk) - sshn   (:,:)
+            gde3w_n(:,:,jk) =   gdept_n(:,:,jk) - ssh    (:,:,Nnn)
          END DO
       ELSE                    ! no ISF cavities 
@@ -375 +376 @@
             gdept_n(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:)
             gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:)
-            gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:)
+            gde3w_n(:,:,jk) = gdept_n(:,:,jk) - ssh    (:,:,Nnn)
          END DO
       ENDIF
@@ -420 +421 @@
 !!gm  Question: use jpdom_data above to read data over jpi x jpj    (like is dom_hgr_read and dom_zgr_read)
 !!              so that it will work with land processor suppression
-!               CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn, ldxios = lrxios    )
-!               CALL iom_get( numror, jpdom_autoglo, 'sshb'   , sshb, ldxios = lrxios    )
+!               CALL iom_get( numror, jpdom_autoglo, 'sshn'   , ssh(:,:,Nnn), ldxios = lrxios    )
+!               CALL iom_get( numror, jpdom_autoglo, 'sshb'   , ssh(:,:,Nbb), ldxios = lrxios    )
 !!gm 
-               CALL iom_get( numror, jpdom_data, 'sshn'   , sshn, ldxios = lrxios    )
-               CALL iom_get( numror, jpdom_data, 'sshb'   , sshb, ldxios = lrxios    )
+               CALL iom_get( numror, jpdom_data, 'sshn'   , ssh(:,:,Nnn), ldxios = lrxios    )
+               CALL iom_get( numror, jpdom_data, 'sshb'   , ssh(:,:,Nbb), ldxios = lrxios    )
 !!gm end
-               IF( l_1st_euler ) THEN
-                  sshb(:,:) = sshn(:,:)
-               ENDIF
+               IF( l_1st_euler )   ssh(:,:,Nbb) = ssh(:,:,Nnn)
             ELSE IF( id1 > 0 ) THEN
                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshn not found in restart files'
-               IF(lwp) write(numout,*) '   set sshn = sshb  and force l_1st_euler = true'
+               IF(lwp) write(numout,*) '   set ssh(Nnn) = ssh(Nbb)  and force l_1st_euler = true'
 !!gm               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
-               CALL iom_get( numror, jpdom_data, 'sshb', sshb, ldxios = lrxios )
-               sshn(:,:) = sshb(:,:)
+               CALL iom_get( numror, jpdom_data, 'sshb', ssh(:,:,Nbb), ldxios = lrxios )
+               ssh(:,:,Nnn) = ssh(:,:,Nbb)
                l_1st_euler = .TRUE.
             ELSE IF( id2 > 0 ) THEN
                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshb not found in restart files'
-               IF(lwp) write(numout,*) 'set sshb = sshn  and force l_1st_euler = true'
-               CALL iom_get( numror, jpdom_data, 'sshn', sshb, ldxios = lrxios )
-               sshb(:,:) = sshn(:,:)
+               IF(lwp) write(numout,*) 'set ssh(Nbb) = ssh(Nnn)  and force l_1st_euler = true'
+               CALL iom_get( numror, jpdom_data, 'sshn', ssh(:,:,Nnn), ldxios = lrxios )
+               ssh(:,:,Nbb) = ssh(:,:,Nnn)
                l_1st_euler = .TRUE.
             ELSE
                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshb and sshn not found in restart file'
-               IF(lwp) write(numout,*) 'set sshb = sshn = 0  and force l_1st_euler = true'
-               sshb(:,:) = 0._wp
-               sshn(:,:) = 0._wp
+               IF(lwp) write(numout,*) 'set ssh(Nbb) = ssh(Nnn) = 0  and force l_1st_euler = true'
+               ssh(:,:,Nbb) = 0._wp
+               ssh(:,:,Nnn) = 0._wp
                l_1st_euler = .TRUE.
             ENDIF
@@ -454 +453 @@
                !
                IF( cn_cfg == 'wad' ) THEN             ! Wetting and drying test case
-                  CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )
-                  tsn  (:,:,:,:) = tsb (:,:,:,:)            ! set now values from to before ones
-                  sshn (:,:)     = sshb(:,:)
-                  un   (:,:,:)   = ub  (:,:,:)
-                  vn   (:,:,:)   = vb  (:,:,:)
+                  CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, ssh(:,:,Nbb)  )
+                  tsn(:,:,:,:) = tsb(:,:,:,:)            ! set now values from to before ones
+                  ssh(:,:,Nnn) = ssh(:,:,Nbb)
+                  un (:,:,:)   = ub (:,:,:)
+                  vn (:,:,:)   = vb (:,:,:)
                ELSE                                   ! Not the test case
-                  sshn(:,:) = -ssh_ref
-                  sshb(:,:) = -ssh_ref
+                  ssh(:,:,Nnn) = -ssh_ref
+                  ssh(:,:,Nbb) = -ssh_ref
                   !
                   DO jj = 1, jpj
                      DO ji = 1, jpi
                         IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN  ! if total depth is less than min depth
-                           sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
-                           sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
-                           ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
+                           ssh(ji,jj,Nbb) = rn_wdmin1 - (ht_0(ji,jj) )
+                           ssh(ji,jj,Nnn) = rn_wdmin1 - (ht_0(ji,jj) )
+                           ssh(ji,jj,Naa) = rn_wdmin1 - (ht_0(ji,jj) )
                         ENDIF
                      END DO
@@ -485 +484 @@
                !
                ! Just to read set ssh in fact, called latter once vertical grid is set up:
-!               CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb  )
-               sshn(:,:) = 0._wp
-               sshb(:,:) = 0._wp
+!               CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, ssh(:,:,Nbb)  )
+               ssh(:,:,Nnn) = 0._wp
+               ssh(:,:,Nbb) = 0._wp
                !
             END IF
@@ -496 +495 @@
          !                                   ! ===================
 
-!!gm      DO NOTHING,   sshb and sshn  are written in restart.F90
+!!gm      DO NOTHING,   ssh(Nbb) and ssh(Nnn)  are written in restart.F90
 
       ENDIF
@@ -594 +593 @@
       DO jj = 1, jpjm1                    ! start from 1 due to f-point
          DO ji = 1, jpim1
-            zsshu_h(ji,jj) = 0.50_wp * ( sshn(ji  ,jj) + sshn(ji+1,jj  ) ) * ssumask(ji,jj)
-            zsshv_h(ji,jj) = 0.50_wp * ( sshn(ji  ,jj) + sshn(ji  ,jj+1) ) * ssvmask(ji,jj)
-            zsshf_h(ji,jj) = 0.25_wp * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)   & 
-               &                       + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) * ssfmask(ji,jj)
+            zsshu_h(ji,jj) = 0.50_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji+1,jj  ,Nnn) ) * ssumask(ji,jj)
+            zsshv_h(ji,jj) = 0.50_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn) ) * ssvmask(ji,jj)
+            zsshf_h(ji,jj) = 0.25_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn)   & 
+               &                       + ssh(ji+1,jj,Nnn) + ssh(ji+1,jj+1,Nnn) ) * ssfmask(ji,jj)
          END DO
       END DO      
@@ -604 +603 @@
       !                             !==  ht, hu and hv  == !   (and their inverse)
       !
-      ht_n   (:,:) = ht_0(:,:) +  sshn  (:,:)
+      ht_n   (:,:) = ht_0(:,:) +  ssh   (:,:,Nnn)
       hu_n   (:,:) = hu_0(:,:) + zsshu_h(:,:)
       hv_n   (:,:) = hv_0(:,:) + zsshv_h(:,:)
@@ -612 +611 @@
       !                             !==  ssh / h  factor at t-, u- ,v- & f-points  ==!
       !
-      zssht_h(:,:) =  sshn  (:,:) * r1_ht_0(:,:)
-      zsshu_h(:,:) = zsshu_h(:,:) * r1_hu_0(:,:)
-      zsshv_h(:,:) = zsshv_h(:,:) * r1_hv_0(:,:)
-      zsshf_h(:,:) = zsshf_h(:,:) * r1_hf_0(:,:)
+      zssht_h(:,:) =  ssh   (:,:,Nnn) * r1_ht_0(:,:)
+      zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)
+      zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)
+      zsshf_h(:,:) = zsshf_h(:,:)     * r1_hf_0(:,:)
       !
       !                             !==  e3t, e3w  ,  e3u, e3uw ,  e3v, e3vw  , and e3f  ==!
       !      
       DO jk = 1, jpkm1
-          e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
-          e3w_n(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
+          e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) *      tmask(:,:,jk)                     )
+          e3w_n(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk) , tmask(:,:,jk+1) ) )
           !
           e3u_n(:,:,jk) =  e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) *  umask(:,:,jk) )
@@ -634 +633 @@
       !                             !== depth of t- and w-points  ==!
       !
-      zssht_h(:,:) = 1._wp + zssht_h(:,:)     ! = 1 + sshn / ht_0
+      zssht_h(:,:) = 1._wp + zssht_h(:,:)     ! = 1 + ssh(Nnn) / ht_0
       !
       IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness 
@@ -640 +639 @@
             gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:)
             gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:)
-            gde3w_n(:,:,jk) =   gdept_n(:,:,jk) - sshn(:,:)
+            gde3w_n(:,:,jk) =   gdept_n(:,:,jk) - ssh(:,:,Nnn)
          END DO
       ELSE                    ! no ISF cavities
+!!gm BUG ???    gdept should be updated down to the ocean floor !  ===>> jpk NOT jpkm1 !!!
          DO jk = 1, jpkm1
             gdept_n(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:)
             gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:)
-            gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:)
+            gde3w_n(:,:,jk) = gdept_n(:,:,jk) - ssh(:,:,Nnn)
          END DO
       ENDIF
@@ -664 +664 @@
       DO jj = 2, jpjm1
          DO ji = 2, jpim1
-            zsshu_h(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji+1,jj  ) ) * ssumask(ji,jj)
-            zsshv_h(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji  ,jj+1) ) * ssvmask(ji,jj)
+            zsshu_h(ji,jj) = 0.5_wp  * ( ssh(ji,jj,Nbb) + ssh(ji+1,jj  ,Nbb) ) * ssumask(ji,jj)
+            zsshv_h(ji,jj) = 0.5_wp  * ( ssh(ji,jj,Nbb) + ssh(ji  ,jj+1,Nbb) ) * ssvmask(ji,jj)
          END DO
       END DO      
@@ -678 +678 @@
       !      
       !                             !==  ssh / h  factor at t-, u- ,v- & f-points  ==!
-      zssht_h(:,:) = sshb (:,:) * r1_ht_0(:,:)
-      zsshu_h   (:,:) = zsshu_h(:,:) * r1_hu_0(:,:)
-      zsshv_h   (:,:) = zsshv_h(:,:) * r1_hv_0(:,:)
+      zssht_h(:,:) = ssh    (:,:,Nbb) * r1_ht_0(:,:)
+      zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)
+      zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)
       !
       !                             !==  e3t, e3w  ,  e3u, e3uw , and  e3v, e3vw  ==!
       DO jk = 1, jpkm1
-          e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
-          e3w_b(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
+          e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) *      tmask(:,:,jk)                     )
+          e3w_b(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk) , tmask(:,:,jk+1) ) )
           !
           e3u_b(:,:,jk) =  e3u_0(:,:,jk) * ( 1._wp + zsshu_h  (:,:) *  umask(:,:,jk) )

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/iscplhsb.F90

@@ -87 +87 @@
       ! diagnose the heat, salt and volume input and compute the correction variable
       !==============================================================================
-
       ! 
-      zdssh(:,:) = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:)
-      IF (.NOT. ln_linssh ) zdssh = 0.0_wp ! already included in the levels by definition
-      
+      zdssh(:,:) = ssh(:,:,Nnn) * ssmask(:,:) - ssh(:,:,Nbb) * psmask_b(:,:)
+      IF (.NOT. ln_linssh )   zdssh = 0._wp     ! already included in the levels by definition
+      !
       DO jk = 1,jpk-1
          DO jj = 2,jpj-1
             DO ji = fs_2,fs_jpim1
-               IF (tmask_h(ji,jj) == 1._wp) THEN
-
-                  ! volume differences
+               IF ( tmask_h(ji,jj) == 1._wp ) THEN
+                  !
+                  !                                   ! volume differences
                   zde3t = e3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk)
 
-                  ! heat diff
-                  zdtem = tsn(ji,jj,jk,jp_tem) * e3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
-                        - tsb(ji,jj,jk,jp_tem) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk)
-                  ! salt diff
-                  zdsal = tsn(ji,jj,jk,jp_sal) * e3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
-                        - tsb(ji,jj,jk,jp_sal) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk)
-               
-                  ! shh changes
+                  !                                   ! heat diff
+                  zdtem = tsn(ji,jj,jk,jp_tem) *  e3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
+                        - tsb(ji,jj,jk,jp_tem) * pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk)
+                  !                                   ! salt diff
+                  zdsal = tsn(ji,jj,jk,jp_sal) *  e3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
+                        - tsb(ji,jj,jk,jp_sal) * pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk)
+                  !                                   ! shh changes
                   IF ( ptmask_b(ji,jj,jk) == 1._wp .OR. tmask(ji,jj,jk) == 1._wp ) THEN 
-                     zde3t = zde3t + zdssh(ji,jj) ! zdssh = 0 if vvl
+                     zde3t = zde3t + zdssh(ji,jj)     ! zdssh = 0 if vvl
                      zdssh(ji,jj) = 0._wp
                   END IF

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/iscplrst.F90

@@ -44 +44 @@
       !! 
       !! ** Purpose : compute initialisation
-      !!              compute extrapolation of restart variable un, vn, tsn, sshn (wetting/drying)   
+      !!              compute extrapolation of restart variable un, vn, tsn, ssh(Nnn) (wetting/drying)   
       !!              compute correction term if needed
       !! 
@@ -91 +91 @@
       !
       !                       ! set _b and _n variables equal
-      tsb (:,:,:,:) = tsn (:,:,:,:)
-      ub  (:,:,:)   = un  (:,:,:)
-      vb  (:,:,:)   = vn  (:,:,:)
-      sshb(:,:)     = sshn(:,:)
+      tsb(:,:,:,:) = tsn(:,:,:,:)
+      ub (:,:,:)   = un (:,:,:)
+      vb (:,:,:)   = vn (:,:,:)
+      ssh(:,:,Nbb) = ssh(:,:,Nnn)
       !
       !                       ! set _b and _n vertical scale factor equal
@@ -117 +117 @@
       !!                   ***  ROUTINE iscpl_rst_interpol  ***
       !! 
-      !! ** Purpose :   compute new tn, sn, un, vn and sshn in case of evolving geometry of ice shelves 
+      !! ** Purpose :   compute new tn, sn, un, vn and ssh(Nnn) in case of evolving geometry of ice shelves 
       !!                compute 2d fields of heat, salt and volume correction
       !! 
@@ -155 +155 @@
       !    
       !                 ! compute new ssh if we open a full water column (average of the closest neigbourgs)  
-      sshb (:,:)=sshn(:,:)
-      zssh0(:,:)=sshn(:,:)
-      zsmask0(:,:) = psmask_b(:,:)
-      zsmask1(:,:) = psmask_b(:,:)
+      ssh  (:,:,Nbb) = ssh(:,:,Nnn)
+      zssh0(:,:)     = ssh(:,:,Nnn)
+      zsmask0(:,:)   = psmask_b(:,:)
+      zsmask1(:,:)   = psmask_b(:,:)
       DO iz = 1, 10                 ! need to be tuned (configuration dependent) (OK for ISOMIP+)
          zdsmask(:,:) = ssmask(:,:)-zsmask0(:,:)
@@ -165 +165 @@
                summsk = zsmask0(ji+1,jj)+zsmask0(ji-1,jj)+zsmask0(ji,jj+1)+zsmask0(ji,jj-1)
                IF (zdsmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
-                  sshn(ji,jj)=( zssh0(ji+1,jj)*zsmask0(ji+1,jj)     &
-                  &           + zssh0(ji-1,jj)*zsmask0(ji-1,jj)     &
-                  &           + zssh0(ji,jj+1)*zsmask0(ji,jj+1)     &
-                  &           + zssh0(ji,jj-1)*zsmask0(ji,jj+1)) / summsk
+                  ssh(ji,jj,Nnn)=( zssh0(ji+1,jj)*zsmask0(ji+1,jj)     &
+                  &              + zssh0(ji-1,jj)*zsmask0(ji-1,jj)     &
+                  &              + zssh0(ji,jj+1)*zsmask0(ji,jj+1)     &
+                  &              + zssh0(ji,jj-1)*zsmask0(ji,jj+1)) / summsk
                   zsmask1(ji,jj) = 1._wp
                ENDIF
             END DO
          END DO
-         CALL lbc_lnk_multi( sshn, 'T', 1., zsmask1, 'T', 1. )
-         zssh0   = sshn
+         CALL lbc_lnk_multi( ssh(:,:,Nnn), 'T', 1., zsmask1, 'T', 1. )
+         zssh0(:,:) = ssh(:,:,Nnn)
          zsmask0 = zsmask1
       END DO
-      sshn(:,:) = sshn(:,:) * ssmask(:,:)
+      ssh(:,:,Nnn) = ssh(:,:,Nnn) * ssmask(:,:)
 
 !=============================================================================
@@ -201 +201 @@
                IF ( tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp ) THEN
                   DO jk = 1, jpk
-                     e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( 1._wp + sshn(ji,jj) * r1_ht_0(ji,jj) * tmask(ji,jj,jk) )
+                     e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( 1._wp + ssh(ji,jj,Nnn) * r1_ht_0(ji,jj) * tmask(ji,jj,jk) )
                   END DO
                ENDIF
@@ -214 +214 @@
          !                                            !* ssh at u- and v-points)
          DO jj = 1, jpjm1   ;   DO ji = 1, jpim1            ! start from 1 due to f-point
-            zsshu(ji,jj) = 0.5_wp  * ( sshn(ji  ,jj) + sshn(ji+1,jj  ) ) * ssumask(ji,jj)
-            zsshv(ji,jj) = 0.5_wp  * ( sshn(ji  ,jj) + sshn(ji  ,jj+1) ) * ssvmask(ji,jj)
-            zsshf(ji,jj) = 0.25_wp * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)   & 
-               &                     + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) * ssfmask(ji,jj)
+            zsshu(ji,jj) = 0.5_wp  * ( ssh(ji  ,jj,Nnn) + ssh(ji+1,jj  ,Nnn) ) * ssumask(ji,jj)
+            zsshv(ji,jj) = 0.5_wp  * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn) ) * ssvmask(ji,jj)
+            zsshf(ji,jj) = 0.25_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn)   & 
+               &                     + ssh(ji+1,jj,Nnn) + ssh(ji+1,jj+1,Nnn) ) * ssfmask(ji,jj)
          END DO             ;   END DO      
          CALL lbc_lnk_multi( zsshu(:,:),'U', 1._wp , zsshu(:,:),'V', 1._wp , zsshf(:,:),'F', 1._wp )
          !
          !                                            !* hu and hv (and their inverse) 
-         ht_n   (:,:) = ht_0(:,:) +  sshn(:,:)
+         ht_n   (:,:) = ht_0(:,:) +  ssh (:,:,Nnn)
          hu_n   (:,:) = hu_0(:,:) + zsshu(:,:)
          hv_n   (:,:) = hv_0(:,:) + zsshv(:,:)
@@ -229 +229 @@
          !
          !                                            !* e3u, e3uw  and  e3v, e3vw
-         z_ssh_h0(:,:) = sshn (:,:) * r1_ht_0(:,:)           ! t-point
+         z_ssh_h0(:,:) = ssh(:,:,Nnn) * r1_ht_0(:,:)         ! t-point
          DO jk = 1, jpkm1
             e3w_n(:,:,jk) = e3w_0(:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * tmask(:,:,jk) )
@@ -248 +248 @@
          END DO
 
-         z_ssh_h0(:,:) = 1._wp + sshn(:,:) * r1_ht_0(:,:)    ! t-point
+         z_ssh_h0(:,:) = 1._wp + ssh(:,:,Nnn) * r1_ht_0(:,:)    ! t-point
          !
          IF( ln_isfcav ) THEN    ! iceshelf cavities : ssh scaling not applied over the iceshelf thickness 
@@ -254 +254 @@
                gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - risfdep(:,:) ) * z_ssh_h0(:,:) + risfdep(:,:)
                gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - risfdep(:,:) ) * z_ssh_h0(:,:) + risfdep(:,:)
-               gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:)
+               gde3w_n(:,:,jk) = gdept_n(:,:,jk) - ssh(:,:,Nnn)
             END DO
          ELSE
@@ -260 +260 @@
                gdept_n(:,:,jk) = gdept_0(:,:,jk) * z_ssh_h0(:,:)
                gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * z_ssh_h0(:,:)
-               gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:)
+               gde3w_n(:,:,jk) = gdept_n(:,:,jk) - ssh(:,:,Nnn)
             END DO
          ENDIF

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/istate.F90

@@ -99 +99 @@
             CALL dta_tsd( nit000, tsb )       ! read 3D T and S data at nit000
             !
-            sshb(:,:)   = 0._wp               ! set the ocean at rest
-            ub  (:,:,:) = 0._wp
-            vb  (:,:,:) = 0._wp  
+            ssh(:,:,Nbb) = 0._wp              ! set the ocean at rest
+            ub (:,:,:)   = 0._wp
+            vb (:,:,:)   = 0._wp  
             !
          ELSE                                 ! user defined initial T and S
-            CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )         
+            CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, ssh(:,:,Nbb) )         
          ENDIF
-         tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones
-         sshn (:,:)     = sshb(:,:)   
-         un   (:,:,:)   = ub  (:,:,:)
-         vn   (:,:,:)   = vb  (:,:,:)
+         tsn  (:,:,:,:) = tsb(:,:,:,:)       ! set now values from to before ones
+         ssh  (:,:,Nnn) = ssh(:,:,Nbb)
+         un   (:,:,:)   = ub (:,:,:)
+         vn   (:,:,:)   = vb (:,:,:)
          hdivn(:,:,jpk) = 0._wp               ! bottom divergence set one for 0 to zero at jpk level
          CALL div_hor( 0 )                    ! compute interior hdivn value  
@@ -115 +115 @@
 
 !!gm POTENTIAL BUG :
-!!gm  ISSUE :  if sshb /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed
+!!gm  ISSUE :  if ssh(Nbb) /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed
 !!             as well as gdept and gdepw....   !!!!! 
 !!      ===>>>>   probably a call to domvvl initialisation here....
@@ -131 +131 @@
          !
 !!gm This is to be changed !!!!
-!         ! - ML - sshn could be modified by istate_eel, so that initialization of e3t_b is done here
+!         ! - ML - ssh(Nnn) could be modified by istate_eel, so that initialization of e3t_b is done here
 !         IF( .NOT.ln_linssh ) THEN
 !            DO jk = 1, jpk

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/restart.F90

@@ -155 +155 @@
          CALL iom_rstput( kt, nitrst, numrow, 'tb'     , tsb(:,:,:,jp_tem), ldxios = lwxios )
          CALL iom_rstput( kt, nitrst, numrow, 'sb'     , tsb(:,:,:,jp_sal), ldxios = lwxios )
-         CALL iom_rstput( kt, nitrst, numrow, 'sshb'   , sshb             , ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'sshb'   , ssh(:,:,Nbb)     , ldxios = lwxios )
          !
          CALL iom_rstput( kt, nitrst, numrow, 'un'     , un               , ldxios = lwxios )     ! now fields
@@ -161 +161 @@
          CALL iom_rstput( kt, nitrst, numrow, 'tn'     , tsn(:,:,:,jp_tem), ldxios = lwxios )
          CALL iom_rstput( kt, nitrst, numrow, 'sn'     , tsn(:,:,:,jp_sal), ldxios = lwxios )
-         CALL iom_rstput( kt, nitrst, numrow, 'sshn'   , sshn             , ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'sshn'   , ssh(:,:,Nnn)     , ldxios = lwxios )
          CALL iom_rstput( kt, nitrst, numrow, 'rhop'   , rhop             , ldxios = lwxios )
          !
@@ -281 +281 @@
       
       IF( iom_varid( numror, 'ub', ldstop = .FALSE. ) > 0  .AND. .NOT.l_1st_euler ) THEN
-         CALL iom_get( numror, jpdom_autoglo, 'ub'     , ub, ldxios = lrxios                )   ! before fields
-         CALL iom_get( numror, jpdom_autoglo, 'vb'     , vb, ldxios = lrxios                )
-         CALL iom_get( numror, jpdom_autoglo, 'tb'     , tsb(:,:,:,jp_tem), ldxios = lrxios )
-         CALL iom_get( numror, jpdom_autoglo, 'sb'     , tsb(:,:,:,jp_sal), ldxios = lrxios )
-         CALL iom_get( numror, jpdom_autoglo, 'sshb'   , sshb, ldxios = lrxios              )
+         CALL iom_get( numror, jpdom_autoglo, 'ub'  , ub               , ldxios = lrxios )   ! before fields
+         CALL iom_get( numror, jpdom_autoglo, 'vb'  , vb               , ldxios = lrxios )
+         CALL iom_get( numror, jpdom_autoglo, 'tb'  , tsb(:,:,:,jp_tem), ldxios = lrxios )
+         CALL iom_get( numror, jpdom_autoglo, 'sb'  , tsb(:,:,:,jp_sal), ldxios = lrxios )
+         CALL iom_get( numror, jpdom_autoglo, 'sshb', ssh(:,:,Nbb)     , ldxios = lrxios )
       ELSE
          l_1st_euler =  .TRUE.      ! before field not found, forced euler 1st time-step
       ENDIF
       !
-      CALL iom_get( numror, jpdom_autoglo, 'un'     , un, ldxios = lrxios )   ! now    fields
-      CALL iom_get( numror, jpdom_autoglo, 'vn'     , vn, ldxios = lrxios )
+      CALL iom_get( numror, jpdom_autoglo, 'un'     , un               , ldxios = lrxios )   ! now    fields
+      CALL iom_get( numror, jpdom_autoglo, 'vn'     , vn               , ldxios = lrxios )
       CALL iom_get( numror, jpdom_autoglo, 'tn'     , tsn(:,:,:,jp_tem), ldxios = lrxios )
       CALL iom_get( numror, jpdom_autoglo, 'sn'     , tsn(:,:,:,jp_sal), ldxios = lrxios )
-      CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn, ldxios = lrxios )
+      CALL iom_get( numror, jpdom_autoglo, 'sshn'   , ssh(:,:,Nnn)     , ldxios = lrxios )
       IF( iom_varid( numror, 'rhop', ldstop = .FALSE. ) > 0 ) THEN
-         CALL iom_get( numror, jpdom_autoglo, 'rhop'   , rhop, ldxios = lrxios )   ! now    potential density
+         CALL iom_get( numror, jpdom_autoglo, 'rhop', rhop             , ldxios = lrxios )   ! now    potential density
       ELSE
          CALL eos( tsn, rhd, rhop, gdept_n(:,:,:) )   
@@ -302 +302 @@
       !
       IF( l_1st_euler ) THEN              ! Euler restart
-         tsb (:,:,:,:) = tsn (:,:,:,:)          ! all before fields set to now values
-         ub  (:,:,:)   = un  (:,:,:)
-         vb  (:,:,:)   = vn  (:,:,:)
-         sshb(:,:)     = sshn(:,:)
+         tsb(:,:,:,:) = tsn(:,:,:,:)          ! all before fields set to now values
+         ub (:,:,:)   = un (:,:,:)
+         vb (:,:,:)   = vn (:,:,:)
+         ssh(:,:,Nbb) = ssh(:,:,Nnn)
          IF( .NOT.ln_linssh )   e3t_b(:,:,:) = e3t_n(:,:,:)
       ENDIF

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dyncor.F90

@@ -1 @@
-MODULE dynvor
+MODULE dyncor
    !!======================================================================
-   !!                       ***  MODULE  dynvor  ***
-   !! Ocean dynamics: Update the momentum trend with the relative and
-   !!                 planetary vorticity trends
+   !!                       ***  MODULE  dyncor  ***
+   !! Ocean dynamics: Update the momentum trend with the planetary vorticity trends
    !!======================================================================
-   !! History :  OPA  ! 1989-12  (P. Andrich)  vor_ens: Original code
-   !!            5.0  ! 1991-11  (G. Madec)  vor_ene, vor_mix: Original code
-   !!            6.0  ! 1996-01  (G. Madec)  s-coord, suppress work arrays
-   !!   NEMO     0.5  ! 2002-08  (G. Madec)  F90: Free form and module
-   !!            1.0  ! 2004-02  (G. Madec)  vor_een: Original code
-   !!             -   ! 2003-08  (G. Madec)  add vor_ctl
-   !!             -   ! 2005-11  (G. Madec)  add dyn_vor (new step architecture)
-   !!            2.0  ! 2006-11  (G. Madec)  flux form advection: add metric term
-   !!            3.2  ! 2009-04  (R. Benshila)  vvl: correction of een scheme
-   !!            3.3  ! 2010-10  (C. Ethe, G. Madec)  reorganisation of initialisation phase
-   !!            3.7  ! 2014-04  (G. Madec)  trend simplification: suppress jpdyn_trd_dat vorticity 
-   !!             -   ! 2014-06  (G. Madec)  suppression of velocity curl from in-core memory
-   !!             -   ! 2016-12  (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T)
-   !!            4.0  ! 2017-07  (G. Madec)  linear dynamics + trends diag. with Stokes-Coriolis
-   !!             -   ! 2018-03  (G. Madec)  add two new schemes (ln_dynvor_enT and ln_dynvor_eet)
-   !!             -   ! 2018-04  (G. Madec)  add pre-computed gradient for metric term calculation
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   dyn_vor       : Update the momentum trend with the vorticity trend
+   !! History :  5.0  ! 2018-07  (G. Madec)  Coriolis trend for Flux Form
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   dyn_cor       : Update the momentum trend with the vorticity trend
    !!       vor_ens   : enstrophy conserving scheme       (ln_dynvor_ens=T)
    !!       vor_ene   : energy conserving scheme          (ln_dynvor_ene=T)
    !!       vor_een   : energy and enstrophy conserving   (ln_dynvor_een=T)
-   !!   dyn_vor_init  : set and control of the different vorticity option
+   !!   dyn_cor_init  : set and control of the different vorticity option
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers
@@ -48 +32 @@
    PRIVATE
 
-   PUBLIC   dyn_vor        ! routine called by step.F90
-   PUBLIC   dyn_vor_init   ! routine called by nemogcm.F90
+   PUBLIC   dyn_cor        ! routine called by step.F90
+   PUBLIC   dyn_cor_init   ! routine called by nemogcm.F90
 
    !                                   !!* Namelist namdyn_vor: vorticity term
@@ -96 +80 @@
 CONTAINS
 
-   SUBROUTINE dyn_vor( kt )
+   SUBROUTINE dyn_cor( kt )
       !!----------------------------------------------------------------------
       !!
@@ -111 +95 @@
       !!----------------------------------------------------------------------
       !
-      IF( ln_timing )   CALL timing_start('dyn_vor')
+      IF( ln_timing )   CALL timing_start('dyn_cor')
       !
       IF( l_trddyn ) THEN     !==  trend diagnostics case : split the added trend in two parts  ==!
@@ -119 +103 @@
          ztrdu(:,:,:) = ua(:,:,:)            !* planetary vorticity trend (including Stokes-Coriolis force)
          ztrdv(:,:,:) = va(:,:,:)
-         SELECT CASE( nvor_scheme )
-         CASE( np_ENS )           ;   CALL vor_ens( kt, ncor, un , vn , ua, va )   ! enstrophy conserving scheme
-            IF( ln_stcor )            CALL vor_ens( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
-         CASE( np_ENE, np_MIX )   ;   CALL vor_ene( kt, ncor, un , vn , ua, va )   ! energy conserving scheme
-            IF( ln_stcor )            CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
-         CASE( np_ENT )           ;   CALL vor_enT( kt, ncor, un , vn , ua, va )   ! energy conserving scheme (T-pts)
-            IF( ln_stcor )            CALL vor_enT( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
-         CASE( np_EET )           ;   CALL vor_eeT( kt, ncor, un , vn , ua, va )   ! energy conserving scheme (een with e3t)
-            IF( ln_stcor )            CALL vor_eeT( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
-         CASE( np_EEN )           ;   CALL vor_een( kt, ncor, un , vn , ua, va )   ! energy & enstrophy scheme
-            IF( ln_stcor )            CALL vor_een( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
-         END SELECT
+
+         CALL cor_ene( kt, ncor, un , vn , ua, va )   ! energy conserving scheme (T-pts)
+
          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
@@ -138 +113 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            SELECT CASE( nvor_scheme )
-            CASE( np_ENT )           ;   CALL vor_enT( kt, nrvm, un , vn , ua, va )  ! energy conserving scheme (T-pts)
-            CASE( np_EET )           ;   CALL vor_eeT( kt, nrvm, un , vn , ua, va )  ! energy conserving scheme (een with e3t)
-            CASE( np_ENE )           ;   CALL vor_ene( kt, nrvm, un , vn , ua, va )  ! energy conserving scheme
-            CASE( np_ENS, np_MIX )   ;   CALL vor_ens( kt, nrvm, un , vn , ua, va )  ! enstrophy conserving scheme
-            CASE( np_EEN )           ;   CALL vor_een( kt, nrvm, un , vn , ua, va )  ! energy & enstrophy scheme
-            END SELECT
+            CALL cor_ene( kt, nrvm, un , vn , ua, va )  ! energy conserving scheme (T-pts)
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
@@ -152 +121 @@
          DEALLOCATE( ztrdu, ztrdv )
          !
-      ELSE              !==  total vorticity trend added to the general trend  ==!
-         !
-         SELECT CASE ( nvor_scheme )      !==  vorticity trend added to the general trend  ==!
-         CASE( np_ENT )                        !* energy conserving scheme  (T-pts)
-                             CALL vor_enT( kt, ntot, un , vn , ua, va )   ! total vorticity trend
-            IF( ln_stcor )   CALL vor_enT( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
-         CASE( np_EET )                        !* energy conserving scheme (een scheme using e3t)
-                             CALL vor_eeT( kt, ntot, un , vn , ua, va )   ! total vorticity trend
-            IF( ln_stcor )   CALL vor_eeT( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
-         CASE( np_ENE )                        !* energy conserving scheme
-                             CALL vor_ene( kt, ntot, un , vn , ua, va )   ! total vorticity trend
-            IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
-         CASE( np_ENS )                        !* enstrophy conserving scheme
-                             CALL vor_ens( kt, ntot, un , vn , ua, va )  ! total vorticity trend
-            IF( ln_stcor )   CALL vor_ens( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend
-         CASE( np_MIX )                        !* mixed ene-ens scheme
-                             CALL vor_ens( kt, nrvm, un , vn , ua, va )   ! relative vorticity or metric trend (ens)
-                             CALL vor_ene( kt, ncor, un , vn , ua, va )   ! planetary vorticity trend (ene)
-            IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
-         CASE( np_EEN )                        !* energy and enstrophy conserving scheme
-                             CALL vor_een( kt, ntot, un , vn , ua, va )   ! total vorticity trend
-            IF( ln_stcor )   CALL vor_een( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
-         END SELECT
+      ELSE              !==  Coriolis (+metric) trend added to the general trend  ==!
+         !
+         !                    !* energy conserving scheme  (T-pts)
+         IF( ln_stcor )   CALL cor_ene( kt, ntot, usd, vsd , ua, va )   !  Stokes drift
+                          CALL cor_ene( kt, ncor, un , vn  , ua, va )   !  
          !
       ENDIF
@@ -182 +133 @@
          &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
-      IF( ln_timing )   CALL timing_stop('dyn_vor')
-      !
-   END SUBROUTINE dyn_vor
-
-
-   SUBROUTINE vor_enT( kt, kvor, pu, pv, pu_rhs, pv_rhs )
+      IF( ln_timing )   CALL timing_stop('dyn_cor')
+      !
+   END SUBROUTINE dyn_cor
+
+
+   SUBROUTINE cor_ene( kt, kvor, pu, pv, pu_rhs, pv_rhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE vor_enT  ***
       !!
-      !! ** Purpose :   Compute the now total vorticity trend and add it to 
-      !!      the general trend of the momentum equation.
+      !! ** Purpose :   Compute the now Coriolis (+ metric term) trend and  
+      !!      add it to the general trend of the momentum equation.
       !!
       !! ** Method  :   Trend evaluated using now fields (centered in time) 
@@ -228 +179 @@
          CASE ( np_COR )                           !* Coriolis (planetary vorticity)
             zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t_n(:,:,jk)
-         CASE ( np_RVO )                           !* relative vorticity
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1
-                  zwz(ji,jj) = (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
-                     &          - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
-               END DO
-            END DO
-            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity 
-               DO jj = 1, jpjm1
-                  DO ji = 1, jpim1
-                     zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
-                  END DO
-               END DO
-            ENDIF
-            CALL lbc_lnk( zwz, 'F', 1. )
-            DO jj = 2, jpj
-               DO ji = 2, jpi   ! vector opt.
-                  zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ) + zwz(ji,jj  )   &
-                     &                  + zwz(ji-1,jj-1) + zwz(ji,jj-1)   ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk)
-               END DO
-            END DO
          CASE ( np_MET )                           !* metric term
             DO jj = 2, jpj
@@ -254 +184 @@
                   zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)   &
                      &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) * e3t_n(ji,jj,jk)
-               END DO
-            END DO
-         CASE ( np_CRV )                           !* Coriolis + relative vorticity
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1                          ! relative vorticity
-                  zwz(ji,jj) = (   e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)   &
-                     &           - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)   ) * r1_e1e2f(ji,jj)
-               END DO
-            END DO
-            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity 
-               DO jj = 1, jpjm1
-                  DO ji = 1, jpim1
-                     zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
-                  END DO
-               END DO
-            ENDIF
-            CALL lbc_lnk( zwz, 'F', 1. )
-            DO jj = 2, jpj
-               DO ji = 2, jpi   ! vector opt.
-                  zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ) + zwz(ji,jj  )    &
-                     &                                 + zwz(ji-1,jj-1) + zwz(ji,jj-1) )  ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk)
                END DO
             END DO
@@ -304 +213 @@
       END DO                                           !   End of slab
       !                                                ! ===============
-   END SUBROUTINE vor_enT
-
-
-   SUBROUTINE vor_ene( kt, kvor, pun, pvn, pua, pva )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE vor_ene  ***
-      !!
-      !! ** Purpose :   Compute the now total vorticity trend and add it to 
-      !!      the general trend of the momentum equation.
-      !!
-      !! ** Method  :   Trend evaluated using now fields (centered in time) 
-      !!       and the Sadourny (1975) flux form formulation : conserves the
-      !!       horizontal kinetic energy.
-      !!         The general trend of momentum is increased due to the vorticity 
-      !!       term which is given by:
-      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f  mi(e1v*e3v vn) ]
-      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f  mj(e2u*e3u un) ]
-      !!       where rvor is the relative vorticity
-      !!
-      !! ** Action : - Update (ua,va) with the now vorticity term trend
-      !!
-      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
-      !!----------------------------------------------------------------------
-      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
-      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend
-      !
-      INTEGER  ::   ji, jj, jk           ! dummy loop indices
-      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz   ! 2D workspace
-      !!----------------------------------------------------------------------
-      !
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-      ENDIF
-      !
-      !                                                ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-         !
-         SELECT CASE( kvor )                 !==  vorticity considered  ==!
-         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
-            zwz(:,:) = ff_f(:,:) 
-         CASE ( np_RVO )                           !* relative vorticity
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
-                     &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
-               END DO
-            END DO
-         CASE ( np_MET )                           !* metric term
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
-                     &       - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
-               END DO
-            END DO
-         CASE ( np_CRV )                           !* Coriolis + relative vorticity
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)      &
-                     &                        - e1u(ji,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
-               END DO
-            END DO
-         CASE ( np_CME )                           !* Coriolis + metric
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
-                     &                     - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
-               END DO
-            END DO
-         CASE DEFAULT                                             ! error
-            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
-         END SELECT
-         !
-         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
-               END DO
-            END DO
-         ENDIF
-
-         IF( ln_sco ) THEN
-            zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk)
-            zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)
-            zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)
-         ELSE
-            zwx(:,:) = e2u(:,:) * pun(:,:,jk)
-            zwy(:,:) = e1v(:,:) * pvn(:,:,jk)
-         ENDIF
-         !                                   !==  compute and add the vorticity term trend  =!
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
-               zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  )
-               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
-               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
-               pua(ji,jj,jk) = pua(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
-               pva(ji,jj,jk) = pva(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
-            END DO  
-         END DO  
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
-   END SUBROUTINE vor_ene
-
-
-   SUBROUTINE vor_ens( kt, kvor, pun, pvn, pua, pva )
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE vor_ens  ***
-      !!
-      !! ** Purpose :   Compute the now total vorticity trend and add it to
-      !!      the general trend of the momentum equation.
-      !!
-      !! ** Method  :   Trend evaluated using now fields (centered in time)
-      !!      and the Sadourny (1975) flux FORM formulation : conserves the
-      !!      potential enstrophy of a horizontally non-divergent flow. the
-      !!      trend of the vorticity term is given by:
-      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ]
-      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f ]  mi-1[ mj(e2u*e3u un) ]
-      !!      Add this trend to the general momentum trend (ua,va):
-      !!          (ua,va) = (ua,va) + ( voru , vorv )
-      !!
-      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
-      !!
-      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
-      !!----------------------------------------------------------------------
-      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
-      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend
-      !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zuav, zvau   ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz, zww   ! 2D workspace
-      !!----------------------------------------------------------------------
-      !
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-      ENDIF
-      !                                                ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-         !
-         SELECT CASE( kvor )                 !==  vorticity considered  ==!
-         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
-            zwz(:,:) = ff_f(:,:) 
-         CASE ( np_RVO )                           !* relative vorticity
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
-                     &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
-               END DO
-            END DO
-         CASE ( np_MET )                           !* metric term
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
-                     &       - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
-               END DO
-            END DO
-         CASE ( np_CRV )                           !* Coriolis + relative vorticity
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)  &
-                     &                        - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
-               END DO
-            END DO
-         CASE ( np_CME )                           !* Coriolis + metric
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
-                     &                     - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
-               END DO
-            END DO
-         CASE DEFAULT                                             ! error
-            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
-         END SELECT
-         !
-         IF( ln_dynvor_msk ) THEN           !==  mask/unmask vorticity ==!
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
-               END DO
-            END DO
-         ENDIF
-         !
-         IF( ln_sco ) THEN                   !==  horizontal fluxes  ==!
-            zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk)
-            zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)
-            zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)
-         ELSE
-            zwx(:,:) = e2u(:,:) * pun(:,:,jk)
-            zwy(:,:) = e1v(:,:) * pvn(:,:,jk)
-         ENDIF
-         !                                   !==  compute and add the vorticity term trend  =!
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zuav = r1_8 * r1_e1u(ji,jj) * (  zwy(ji  ,jj-1) + zwy(ji+1,jj-1)  &
-                  &                           + zwy(ji  ,jj  ) + zwy(ji+1,jj  )  )
-               zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  &
-                  &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1)  )
-               pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
-               pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
-            END DO  
-         END DO  
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
-   END SUBROUTINE vor_ens
-
-
-   SUBROUTINE vor_een( kt, kvor, pun, pvn, pua, pva )
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE vor_een  ***
-      !!
-      !! ** Purpose :   Compute the now total vorticity trend and add it to 
-      !!      the general trend of the momentum equation.
-      !!
-      !! ** Method  :   Trend evaluated using now fields (centered in time) 
-      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves 
-      !!      both the horizontal kinetic energy and the potential enstrophy
-      !!      when horizontal divergence is zero (see the NEMO documentation)
-      !!      Add this trend to the general momentum trend (ua,va).
-      !!
-      !! ** Action : - Update (ua,va) with the now vorticity term trend
-      !!
-      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
-      !!----------------------------------------------------------------------
-      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
-      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend
-      !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      INTEGER  ::   ierr         ! local integer
-      REAL(wp) ::   zua, zva     ! local scalars
-      REAL(wp) ::   zmsk, ze3f   ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj) ::   zwx , zwy , zwz , z1_e3f
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztnw, ztne, ztsw, ztse
-      !!----------------------------------------------------------------------
-      !
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-      ENDIF
-      !
-      !                                                ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-         !
-         SELECT CASE( nn_een_e3f )           ! == reciprocal of e3 at F-point
-         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  ze3f = (  e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
-                     &    + e3t_n(ji,jj  ,jk)*tmask(ji,jj  ,jk) + e3t_n(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)  )
-                  IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4._wp / ze3f
-                  ELSE                       ;   z1_e3f(ji,jj) = 0._wp
-                  ENDIF
-               END DO
-            END DO
-         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  ze3f = (  e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
-                     &    + e3t_n(ji,jj  ,jk)*tmask(ji,jj  ,jk) + e3t_n(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)  )
-                  zmsk = (                    tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
-                     &                      + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk)  )
-                  IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = zmsk / ze3f
-                  ELSE                       ;   z1_e3f(ji,jj) = 0._wp
-                  ENDIF
-               END DO
-            END DO
-         END SELECT
-         !
-         SELECT CASE( kvor )                 !==  vorticity considered  ==!
-         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ff_f(ji,jj) * z1_e3f(ji,jj)
-               END DO
-            END DO
-         CASE ( np_RVO )                           !* relative vorticity
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ( e2v(ji+1,jj  ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)  &
-                     &         - e1u(ji  ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)
-               END DO
-            END DO
-         CASE ( np_MET )                           !* metric term
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = (   ( pvn(ji+1,jj,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
-                     &           - ( pun(ji,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
-               END DO
-            END DO
-         CASE ( np_CRV )                           !* Coriolis + relative vorticity
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)      &
-                     &                           - e1u(ji  ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  )   &
-                     &                        * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
-               END DO
-            END DO
-         CASE ( np_CME )                           !* Coriolis + metric
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = (   ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
-                     &                         - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
-               END DO
-            END DO
-         CASE DEFAULT                                             ! error
-            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
-         END SELECT
-         !
-         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
-               END DO
-            END DO
-         ENDIF
-         !
-         CALL lbc_lnk( zwz, 'F', 1. )
-         !
-         !                                   !==  horizontal fluxes  ==!
-         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)
-         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)
-
-         !                                   !==  compute and add the vorticity term trend  =!
-         jj = 2
-         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
-         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
-               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
-               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
-               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
-               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
-         END DO
-         DO jj = 3, jpj
-            DO ji = fs_2, jpi   ! vector opt. ok because we start at jj = 3
-               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
-               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
-               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
-               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
-            END DO
-         END DO
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
-                  &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
-               zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
-                  &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
-               pua(ji,jj,jk) = pua(ji,jj,jk) + zua
-               pva(ji,jj,jk) = pva(ji,jj,jk) + zva
-            END DO  
-         END DO  
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
-   END SUBROUTINE vor_een
-
-
-
-   SUBROUTINE vor_eeT( kt, kvor, pun, pvn, pua, pva )
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE vor_eeT  ***
-      !!
-      !! ** Purpose :   Compute the now total vorticity trend and add it to 
-      !!      the general trend of the momentum equation.
-      !!
-      !! ** Method  :   Trend evaluated using now fields (centered in time) 
-      !!      and the Arakawa and Lamb (1980) vector form formulation using 
-      !!      a modified version of Arakawa and Lamb (1980) scheme (see vor_een).
-      !!      The change consists in 
-      !!      Add this trend to the general momentum trend (ua,va).
-      !!
-      !! ** Action : - Update (ua,va) with the now vorticity term trend
-      !!
-      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
-      !!----------------------------------------------------------------------
-      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
-      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend
-      !
-      INTEGER  ::   ji, jj, jk     ! dummy loop indices
-      INTEGER  ::   ierr           ! local integer
-      REAL(wp) ::   zua, zva       ! local scalars
-      REAL(wp) ::   zmsk, z1_e3t   ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj) ::   zwx , zwy , zwz
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztnw, ztne, ztsw, ztse
-      !!----------------------------------------------------------------------
-      !
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-      ENDIF
-      !
-      !                                                ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-         !
-         !
-         SELECT CASE( kvor )                 !==  vorticity considered  ==!
-         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ff_f(ji,jj)
-               END DO
-            END DO
-         CASE ( np_RVO )                           !* relative vorticity
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
-                     &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) &
-                     &       * r1_e1e2f(ji,jj)
-               END DO
-            END DO
-         CASE ( np_MET )                           !* metric term
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
-                     &       - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
-               END DO
-            END DO
-         CASE ( np_CRV )                           !* Coriolis + relative vorticity
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
-                     &                           - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) &
-                     &                      * r1_e1e2f(ji,jj)    )
-               END DO
-            END DO
-         CASE ( np_CME )                           !* Coriolis + metric
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
-                     &                     - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
-               END DO
-            END DO
-         CASE DEFAULT                                             ! error
-            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
-         END SELECT
-         !
-         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
-               END DO
-            END DO
-         ENDIF
-         !
-         CALL lbc_lnk( zwz, 'F', 1. )
-         !
-         !                                   !==  horizontal fluxes  ==!
-         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)
-         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)
-
-         !                                   !==  compute and add the vorticity term trend  =!
-         jj = 2
-         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
-         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
-               z1_e3t = 1._wp / e3t_n(ji,jj,jk)
-               ztne(ji,jj) = ( zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1) ) * z1_e3t
-               ztnw(ji,jj) = ( zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  ) ) * z1_e3t
-               ztse(ji,jj) = ( zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1) ) * z1_e3t
-               ztsw(ji,jj) = ( zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  ) ) * z1_e3t
-         END DO
-         DO jj = 3, jpj
-            DO ji = fs_2, jpi   ! vector opt. ok because we start at jj = 3
-               z1_e3t = 1._wp / e3t_n(ji,jj,jk)
-               ztne(ji,jj) = ( zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1) ) * z1_e3t
-               ztnw(ji,jj) = ( zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  ) ) * z1_e3t
-               ztse(ji,jj) = ( zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1) ) * z1_e3t
-               ztsw(ji,jj) = ( zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  ) ) * z1_e3t
-            END DO
-         END DO
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
-                  &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
-               zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
-                  &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
-               pua(ji,jj,jk) = pua(ji,jj,jk) + zua
-               pva(ji,jj,jk) = pva(ji,jj,jk) + zva
-            END DO  
-         END DO  
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
-   END SUBROUTINE vor_eeT
-
-
-   SUBROUTINE dyn_vor_init
+   END SUBROUTINE cor_ene
+
+
+   SUBROUTINE dyn_cor_init
       !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE dyn_vor_init  ***
+      !!                  ***  ROUTINE dyn_cor_init  ***
       !!
       !! ** Purpose :   Control the consistency between cpp options for
@@ -882 +292 @@
       SELECT CASE( n_dynadv )
       CASE( np_LIN_dyn )
-         IF(lwp) WRITE(numout,*) '   ==>>>   linear dynamics : total vorticity = Coriolis'
+         IF(lwp) WRITE(numout,*) '   ==>>>   linear dynamics : only Coriolis, no metric term'
          nrvm = np_COR        ! planetary vorticity
          ntot = np_COR        !     -         -
       CASE( np_VEC_c2  )
-         IF(lwp) WRITE(numout,*) '   ==>>>   vector form dynamics : total vorticity = Coriolis + relative vorticity' 
-         nrvm = np_RVO        ! relative vorticity
-         ntot = np_CRV        ! relative + planetary vorticity         
+         CALL ctl_stop( 'dyncor_init : cor_ene requires FLUX form dynamics, not VECTOR form' )
       CASE( np_FLX_c2 , np_FLX_ubs  )
          IF(lwp) WRITE(numout,*) '   ==>>>   flux form dynamics : total vorticity = Coriolis + metric term'
@@ -894 +302 @@
          ntot = np_CME        ! Coriolis + metric term
          !
-         SELECT CASE( nvor_scheme )    ! pre-computed gradients for the metric term:
-         CASE( np_ENT )                      !* T-point metric term :   pre-compute di(e2u)/2 and dj(e1v)/2
+         !                    ! pre-computed gradients for the metric term:
+         !                         !* T-point metric term :   pre-compute di(e2u)/2 and dj(e1v)/2
             ALLOCATE( di_e2u_2(jpi,jpj), dj_e1v_2(jpi,jpj) )
             DO jj = 2, jpjm1
@@ -905 +313 @@
             CALL lbc_lnk_multi( di_e2u_2, 'T', -1. , dj_e1v_2, 'T', -1. )   ! Lateral boundary conditions
             !
-         CASE DEFAULT                        !* F-point metric term :   pre-compute di(e2u)/(2*e1e2f) and dj(e1v)/(2*e1e2f)
-            ALLOCATE( di_e2v_2e1e2f(jpi,jpj), dj_e1u_2e1e2f(jpi,jpj) )
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1
-                  di_e2v_2e1e2f(ji,jj) = ( e2v(ji+1,jj  ) - e2v(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
-                  dj_e1u_2e1e2f(ji,jj) = ( e1u(ji  ,jj+1) - e1u(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
-               END DO
-            END DO
-            CALL lbc_lnk_multi( di_e2v_2e1e2f, 'F', -1. , dj_e1u_2e1e2f, 'F', -1. )   ! Lateral boundary conditions
-         END SELECT
          !
       END SELECT
-      
-      IF(lwp) THEN                   ! Print the choice
-         WRITE(numout,*)
-         SELECT CASE( nvor_scheme )
-         CASE( np_ENS )   ;   WRITE(numout,*) '   ==>>>   enstrophy conserving scheme (ENS)'
-         CASE( np_ENE )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at F-points) (ENE)'
-         CASE( np_ENT )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at T-points) (ENT)'
-         CASE( np_EET )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3t) (EET)'
-         CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)'
-         CASE( np_MIX )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)'
-         END SELECT         
-      ENDIF
-      !
-   END SUBROUTINE dyn_vor_init
+      !      
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at F-points) (ENE)'
+      !
+   END SUBROUTINE dyn_cor_init
 
    !!==============================================================================
-END MODULE dynvor
+END MODULE dyncor

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynhpg.F90

@@ -453 +453 @@
         DO jj = 2, jpjm1
            DO ji = 2, jpim1 
-             ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji+1,jj) ) >                &
-                  &    MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            &
-                  &    MAX(  sshn(ji,jj) +  ht_0(ji,jj),  sshn(ji+1,jj) + ht_0(ji+1,jj) )  &
-                  &                                                       > rn_wdmin1 + rn_wdmin2
-             ll_tmp2 = ( ABS( sshn(ji,jj)              -  sshn(ji+1,jj) ) > 1.E-12 ) .AND. (       &
-                  &    MAX(   sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                &
-                  &    MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
+             ll_tmp1 = MIN(  ssh(ji,jj,Nnn)               , ssh(ji+1,jj,Nnn)                  )                            &
+                  &  > MAX(                 - ht_0(ji,jj) ,                  - ht_0(ji+1,jj)  )  .AND.                     &
+                  &    MAX(  ssh(ji,jj,Nnn) + ht_0(ji,jj) , ssh(ji+1,jj,Nnn) + ht_0(ji+1,jj)  )    >   rn_wdmin1 + rn_wdmin2
+                  !
+             ll_tmp2 = ABS( ssh(ji,jj,Nnn) -  ssh(ji+1,jj,Nnn) ) > 1.E-12                        .AND.                     &
+                  &    MAX( ssh(ji,jj,Nnn) ,  ssh(ji+1,jj,Nnn) ) > MAX( ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2
 
              IF(ll_tmp1) THEN
                zcpx(ji,jj) = 1.0_wp
              ELSE IF(ll_tmp2) THEN
-               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
-               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &
-                           &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) )
+               ! no worries about  ssh(ji+1,jj,,Nnn) -  ssh(ji,jj,Nnn) = 0, it won't happen ! here
+               zcpx(ji,jj) = ABS(   ( ssh(ji+1,jj,Nnn) + ht_0(ji+1,jj) - ssh(ji,jj,Nnn) - ht_0(ji,jj) )  &
+                           &      / ( ssh(ji+1,jj,Nnn)                 - ssh(ji,jj,Nnn)               )  )
              ELSE
                zcpx(ji,jj) = 0._wp
              END IF
       
-             ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                &
-                  &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            &
-                  &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji,jj+1) + ht_0(ji,jj+1) )  &
-                  &                                                      > rn_wdmin1 + rn_wdmin2
-             ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji,jj+1) ) > 1.E-12 ) .AND. (        &
-                  &    MAX(   sshn(ji,jj)             ,  sshn(ji,jj+1) ) >                &
-                  &    MAX(  -ht_0(ji,jj)             , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
+             ll_tmp1 = MIN(  ssh (ji,jj,Nnn)             ,  ssh (ji,jj+1,Nnn) )                     &
+                  &  > MAX( -ht_0(ji,jj)                 , -ht_0(ji,jj+1) )              .AND.      &
+                  &    MAX(  ssh(ji,jj,Nnn) + ht_0(ji,jj),  ssh (ji,jj+1,Nnn) + ht_0(ji,jj+1) )     &
+                  &  > rn_wdmin1 + rn_wdmin2
+             ll_tmp2 = ( ABS( ssh (ji,jj,Nnn)            -  ssh (ji,jj+1,Nnn) ) > 1.E-12 ) .AND.    &
+                  &  ( MAX(   ssh (ji,jj,Nnn)            ,  ssh (ji,jj+1,Nnn) ) >                   &
+                  &    MAX(  -ht_0(ji,jj)                , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
 
              IF(ll_tmp1) THEN
                zcpy(ji,jj) = 1.0_wp
              ELSE IF(ll_tmp2) THEN
-               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
-               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &
-                           &    / (sshn(ji,jj+1) - sshn(ji,jj  )) )
+               ! no worries about  ssh(ji,jj+1,Nnn) - ssh(ji,jj  ,Nnn) = 0, it won't happen ! here
+               zcpy(ji,jj) = ABS( ( ssh(ji,jj+1,Nnn) + ht_0(ji,jj+1) - ssh(ji,jj,Nnn) - ht_0(ji,jj) )  &
+                  &             / ( ssh(ji,jj+1,Nnn)                 - ssh(ji,jj,Nnn)               )  )
              ELSE
                zcpy(ji,jj) = 0._wp
@@ -687 +686 @@
         DO jj = 2, jpjm1
            DO ji = 2, jpim1 
-             ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                &
-                  &    MAX( -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) .AND.            &
-                  &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji+1,jj) + ht_0(ji+1,jj) )  &
-                  &                                                      > rn_wdmin1 + rn_wdmin2
-             ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji+1,jj) ) > 1.E-12 ) .AND. (        &
-                  &    MAX(   sshn(ji,jj)             ,  sshn(ji+1,jj) ) >                &
+             ll_tmp1 = MIN(  ssh (ji,jj,Nnn)              ,  ssh (ji+1,jj,Nnn) )                  &
+                  &  > MAX( -ht_0(ji,jj)                  , -ht_0(ji+1,jj)     ) .AND.            &
+                  &    MAX(  ssh (ji,jj,Nnn) + ht_0(ji,jj),  ssh (ji+1,jj,Nnn) + ht_0(ji+1,jj) )  &
+                  & > rn_wdmin1 + rn_wdmin2
+             ll_tmp2 = ( ABS( ssh (ji,jj,Nnn)             -  ssh (ji+1,jj,Nnn) ) > 1.E-12 ) .AND. (   &
+                  &    MAX(   ssh (ji,jj,Nnn)             ,  ssh (ji+1,jj,Nnn) ) >                &
                   &    MAX(  -ht_0(ji,jj)             , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
              IF(ll_tmp1) THEN
                zcpx(ji,jj) = 1.0_wp
              ELSE IF(ll_tmp2) THEN
-               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
-               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &
-                           &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) )
+               ! no worries about  ssh(ji+1,jj,Nnn) -  ssh(ji  ,jj,Nnn) = 0, it won't happen ! here
+               zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Nnn) + ht_0(ji+1,jj) - ssh(ji,jj,Nnn) - ht_0(ji,jj) )  &
+                           &    / (ssh(ji+1,jj,Nnn)                 - ssh(ji,jj,Nnn)               )  )
              ELSE
                zcpx(ji,jj) = 0._wp
              END IF
       
-             ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                &
-                  &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            &
-                  &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji,jj+1) + ht_0(ji,jj+1) )  &
-                  &                                                      > rn_wdmin1 + rn_wdmin2
-             ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji,jj+1) ) > 1.E-12 ) .AND. (        &
-                  &    MAX(   sshn(ji,jj)             ,  sshn(ji,jj+1) ) >                &
+             ll_tmp1 = MIN(  ssh (ji,jj,Nnn)              ,  ssh (ji,jj+1,Nnn) ) >                &
+                  &    MAX( -ht_0(ji,jj)                  , -ht_0(ji,jj+1) )          .AND.       &
+                  &    MAX(  ssh (ji,jj,Nnn) + ht_0(ji,jj),  ssh (ji,jj+1,Nnn) + ht_0(ji,jj+1) )  &
+                  &  > rn_wdmin1 + rn_wdmin2
+             ll_tmp2 = ( ABS( ssh (ji,jj,Nnn)             -  ssh (ji,jj+1,Nnn) ) > 1.E-12 ) .AND. (        &
+                  &    MAX(   ssh (ji,jj,Nnn)             ,  ssh (ji,jj+1,Nnn) ) >                &
                   &    MAX(  -ht_0(ji,jj)             , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
 
@@ -715 +714 @@
                zcpy(ji,jj) = 1.0_wp
              ELSE IF(ll_tmp2) THEN
-               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
-               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &
-                           &    / (sshn(ji,jj+1) - sshn(ji,jj  )) )
+               ! no worries about  ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) = 0, it won't happen ! here
+               zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Nnn) + ht_0(ji,jj+1) - ssh(ji,jj,Nnn) - ht_0(ji,jj) )   &
+                           &    / (ssh(ji,jj+1,Nnn)                 - ssh(ji,jj,Nnn)               )   )
              ELSE
                zcpy(ji,jj) = 0._wp
@@ -975 +974 @@
          DO jj = 2, jpjm1
            DO ji = 2, jpim1 
-             ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                &
+             ll_tmp1 = MIN(  ssh (ji,jj,Nnn)              ,  ssh(ji+1,jj,Nnn) ) >                &
                   &    MAX( -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) .AND.            &
-                  &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji+1,jj) + ht_0(ji+1,jj) )  &
-                  &                                                      > rn_wdmin1 + rn_wdmin2
-             ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji+1,jj) ) > 1.E-12 ) .AND. (         &
-                  &    MAX(   sshn(ji,jj)             ,  sshn(ji+1,jj) ) >                &
+                  &    MAX(  ssh(ji,jj,Nnn) + ht_0(ji,jj),  ssh(ji+1,jj,Nnn) + ht_0(ji+1,jj) )  &
+                  &  > rn_wdmin1 + rn_wdmin2
+             ll_tmp2 = ( ABS( ssh(ji,jj,Nnn)             -  ssh (ji+1,jj,Nnn) ) > 1.E-12 ) .AND. (         &
+                  &    MAX(   ssh(ji,jj,Nnn)             ,  ssh (ji+1,jj,Nnn) ) >                &
                   &    MAX(  -ht_0(ji,jj)             , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
 
@@ -986 +985 @@
                zcpx(ji,jj) = 1.0_wp
              ELSE IF(ll_tmp2) THEN
-               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
-               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &
-                           &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) )
+               ! no worries about  ssh(ji+1,jj,Nnn) -  ssh(ji,jj,Nnn) = 0, it won't happen ! here
+               zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Nnn) + ht_0(ji+1,jj) - ssh(ji,jj,Nnn) - ht_0(ji,jj) )  &
+                           &    / (ssh(ji+1,jj,Nnn)                 - ssh(ji,jj,Nnn)               )  )
               
                 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
@@ -995 +994 @@
              END IF
       
-             ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                &
+             ll_tmp1 = MIN(  ssh (ji,jj,Nnn)          ,  ssh(ji,jj+1,Nnn) ) >                &
                   &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            &
-                  &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji,jj+1) + ht_0(ji,jj+1) )  &
-                  &                                                      > rn_wdmin1 + rn_wdmin2
-             ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji,jj+1) ) > 1.E-12 ) .AND. (      &
-                  &    MAX(   sshn(ji,jj)             ,  sshn(ji,jj+1) ) >                &
+                  &    MAX(  ssh (ji,jj,Nnn) + ht_0(ji,jj),  ssh(ji,jj+1,Nnn) + ht_0(ji,jj+1) )  &
+                  &  > rn_wdmin1 + rn_wdmin2
+             ll_tmp2 = ( ABS( ssh (ji,jj,Nnn)         -  ssh (ji,jj+1,Nnn) ) > 1.E-12 ) .AND. (      &
+                  &    MAX(   ssh (ji,jj,Nnn)         ,  ssh (ji,jj+1,Nnn) ) >                &
                   &    MAX(  -ht_0(ji,jj)             , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
 
@@ -1006 +1005 @@
                zcpy(ji,jj) = 1.0_wp
              ELSE IF(ll_tmp2) THEN
-               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
-               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &
-                           &    / (sshn(ji,jj+1) - sshn(ji,jj  )) )
+               ! no worries about  ssh(ji,jj+1,Nnn) -  ssh(ji,jj,Nnn) = 0, it won't happen ! here
+               zcpy(ji,jj) = ABS(  ( ssh(ji,jj+1,Nnn) + ht_0(ji,jj+1) - ssh(ji,jj,Nnn) - ht_0(ji,jj))   &
+                           &     / ( ssh(ji,jj+1,Nnn)                 - ssh(ji,jj,Nnn)              )   )
                 zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp)
 
@@ -1041 +1040 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            zdept(ji,jj,1) = 0.5_wp * e3w_n(ji,jj,1) - sshn(ji,jj) * znad
+            zdept(ji,jj,1) = 0.5_wp * e3w_n(ji,jj,1) - ssh(ji,jj,Nnn) * znad
          END DO
       END DO
@@ -1087 +1086 @@
 
       ! Prepare zsshu_n and zsshv_n
+!!gm Vector form   
+!         DO jj = 2, jpjm1
+!           DO ji = 2, jpim1
+!   !!gm BUG ?    if it is ssh at u- & v-point then it should be:
+!   !          zsshu_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Nnn) + e1e2t(ji+1,jj) * ssh(ji+1,jj,Nnn)) * &
+!   !                         & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 
+!   !          zsshv_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Nnn) + e1e2t(ji,jj+1) * ssh(ji,jj+1,Nnn)) * &
+!   !                         & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 
+!   !!gm not this:
+!             zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh(ji,jj,Nnn) + e1e2u(ji+1, jj) * ssh(ji+1,jj,Nnn)) * &
+!                            & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 
+!             zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh(ji,jj,Nnn) + e1e2v(ji+1, jj) * ssh(ji,jj+1,Nnn)) * &
+!                            & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 
+!           END DO
+!         END DO
+!!gm Flux form :
       DO jj = 2, jpjm1
-        DO ji = 2, jpim1
-!!gm BUG ?    if it is ssh at u- & v-point then it should be:
-!          zsshu_n(ji,jj) = (e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji+1,jj) * sshn(ji+1,jj)) * &
-!                         & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 
-!          zsshv_n(ji,jj) = (e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji,jj+1) * sshn(ji,jj+1)) * &
-!                         & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 
-!!gm not this:
-          zsshu_n(ji,jj) = (e1e2u(ji,jj) * sshn(ji,jj) + e1e2u(ji+1, jj) * sshn(ji+1,jj)) * &
-                         & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 
-          zsshv_n(ji,jj) = (e1e2v(ji,jj) * sshn(ji,jj) + e1e2v(ji+1, jj) * sshn(ji,jj+1)) * &
-                         & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 
+         DO ji = 2, jpim1
+            zsshu_n(ji,jj) = 0.5_wp * ( ssh(ji,jj,Nnn) + ssh(ji+1,jj,Nnn) ) * ssumask(ji,jj)
+            zsshv_n(ji,jj) = 0.5_wp * ( ssh(ji,jj,Nnn) + ssh(ji,jj+1,Nnn) ) * ssvmask(ji,jj)
         END DO
       END DO
-
       CALL lbc_lnk_multi (zsshu_n, 'U', 1., zsshv_n, 'V', 1. )
-
+      !
       DO jj = 2, jpjm1
-        DO ji = 2, jpim1
-          zu(ji,jj,1) = - ( e3u_n(ji,jj,1) - zsshu_n(ji,jj) * znad) 
-          zv(ji,jj,1) = - ( e3v_n(ji,jj,1) - zsshv_n(ji,jj) * znad)
+         DO ji = 2, jpim1
+            zu(ji,jj,1) = - ( e3u_n(ji,jj,1) - zsshu_n(ji,jj) * znad) 
+            zv(ji,jj,1) = - ( e3v_n(ji,jj,1) - zsshv_n(ji,jj) * znad)
         END DO
       END DO
-
+      !
       DO jk = 2, jpkm1
-        DO jj = 2, jpjm1
-          DO ji = 2, jpim1
-            zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u_n(ji,jj,jk)
-            zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v_n(ji,jj,jk)
-          END DO
-        END DO
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u_n(ji,jj,jk)
+               zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v_n(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u_n(ji,jj,jk)
+               zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v_n(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               zu(ji,jj,jk) = MIN(  zu(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) )  )
+               zu(ji,jj,jk) = MAX(  zu(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) )  )
+               zv(ji,jj,jk) = MIN(  zv(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) )  )
+               zv(ji,jj,jk) = MAX(  zv(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) )  )
+            END DO
+         END DO
       END DO
 
@@ -1123 +1149 @@
         DO jj = 2, jpjm1
           DO ji = 2, jpim1
-            zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u_n(ji,jj,jk)
-            zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v_n(ji,jj,jk)
-          END DO
-        END DO
-      END DO
-
-      DO jk = 1, jpkm1
-        DO jj = 2, jpjm1
-          DO ji = 2, jpim1
-            zu(ji,jj,jk) = MIN(  zu(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) )  )
-            zu(ji,jj,jk) = MAX(  zu(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) )  )
-            zv(ji,jj,jk) = MIN(  zv(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) )  )
-            zv(ji,jj,jk) = MAX(  zv(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) )  )
-          END DO
-        END DO
-      END DO
-
-
-      DO jk = 1, jpkm1
-        DO jj = 2, jpjm1
-          DO ji = 2, jpim1
-            zpwes = 0._wp; zpwed = 0._wp
-            zpnss = 0._wp; zpnsd = 0._wp
+            zpwes = 0._wp   ;   zpwed = 0._wp
+            zpnss = 0._wp   ;   zpnsd = 0._wp
             zuijk = zu(ji,jj,jk)
             zvijk = zv(ji,jj,jk)
@@ -1176 +1181 @@
                DO WHILE ( -zdept(jid,jj,jk1) < zuijk )
                  IF( jk1 == 1 ) THEN
-                   zdeps = zdept(jid,jj,1) + MIN(zuijk, sshn(jid,jj)*znad)
+                   zdeps = zdept(jid,jj,1) + MIN(zuijk, ssh(jid,jj,Nnn)*znad)
                    zrhdt1 = zrhh(jid,jj,1) - interp3(zdept(jid,jj,1), asp(jid,jj,1), &
                                                      bsp(jid,jj,1),   csp(jid,jj,1), &
@@ -1196 +1201 @@
                IF( .NOT.ln_linssh ) THEN
                  zdpdx2 = zcoef0 * r1_e1u(ji,jj) * &
-                    &    ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) )
+                    &    ( REAL(jis-jid, wp) * (zpwes + zpwed) + (ssh(ji+1,jj,Nnn)-ssh(ji,jj,Nnn)) )
                 ELSE
                  zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed)
@@ -1234 +1239 @@
                DO WHILE ( -zdept(ji,jjd,jk1) < zvijk )
                  IF( jk1 == 1 ) THEN
-                   zdeps = zdept(ji,jjd,1) + MIN(zvijk, sshn(ji,jjd)*znad)
+                   zdeps = zdept(ji,jjd,1) + MIN(zvijk, ssh(ji,jjd,Nnn)*znad)
                    zrhdt1 = zrhh(ji,jjd,1) - interp3(zdept(ji,jjd,1), asp(ji,jjd,1), &
                                                      bsp(ji,jjd,1),   csp(ji,jjd,1), &
@@ -1255 +1260 @@
                IF( .NOT.ln_linssh ) THEN
                   zdpdy2 = zcoef0 * r1_e2v(ji,jj) * &
-                          ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) )
+                          ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (ssh(ji,jj+1,Nnn)-ssh(ji,jj,Nnn)) )
                ELSE
                   zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd )

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynnxt.F90

@@ -252 +252 @@
                !                                            !* ssh at u- and v-points)
                DO jj = 2, jpjm1   ;   DO ji = 2, jpim1
-                  zsshu(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji+1,jj  ) ) * ssumask(ji,jj)
-                  zsshv(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji  ,jj+1) ) * ssvmask(ji,jj)
+                  zsshu(ji,jj) = 0.5_wp * ( ssh(ji,jj,Nbb) + ssh(ji+1,jj  ,Nbb) ) * ssumask(ji,jj)
+                  zsshv(ji,jj) = 0.5_wp * ( ssh(ji,jj,Nbb) + ssh(ji  ,jj+1,Nbb) ) * ssvmask(ji,jj)
                END DO             ;   END DO      
                CALL lbc_lnk_multi( zsshu(:,:),'U', 1._wp , zsshu(:,:),'V', 1._wp )
@@ -288 +288 @@
                !                                            !* ssh at u- and v-points)
                DO jj = 2, jpjm1   ;   DO ji = 2, jpim1
-                  zsshu(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji+1,jj  ) ) * ssumask(ji,jj)
-                  zsshv(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji  ,jj+1) ) * ssvmask(ji,jj)
+                  zsshu(ji,jj) = 0.5_wp * ( ssh(ji,jj,Nbb) + ssh(ji+1,jj  ,Nbb) ) * ssumask(ji,jj)
+                  zsshv(ji,jj) = 0.5_wp * ( ssh(ji,jj,Nbb) + ssh(ji  ,jj+1,Nbb) ) * ssvmask(ji,jj)
                END DO             ;   END DO      
                CALL lbc_lnk_multi( zsshu(:,:),'U', 1._wp , zsshu(:,:),'V', 1._wp )
@@ -355 +355 @@
       IF(.NOT.ln_linssh ) THEN
          DO jj = 2, jpjm1   ;   DO ji = 2, jpim1
-            zsshu(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji+1,jj  ) ) * ssumask(ji,jj)
-            zsshv(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji  ,jj+1) ) * ssvmask(ji,jj)
+            zsshu(ji,jj) = 0.5_wp * ( ssh(ji,jj,Nbb) + ssh(ji+1,jj  ,Nbb) ) * ssumask(ji,jj)
+            zsshv(ji,jj) = 0.5_wp * ( ssh(ji,jj,Nbb) + ssh(ji  ,jj+1,Nbb) ) * ssvmask(ji,jj)
          END DO             ;   END DO      
          CALL lbc_lnk_multi( zsshu(:,:),'U', 1._wp , zsshu(:,:),'V', 1._wp )

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynspg.F90

@@ -120 +120 @@
                   DO jj = 2, jpjm1                    ! add tide potential + scalar approximation of load potential
                      DO ji = fs_2, fs_jpim1   ! vector opt.
-                        spgu(ji,jj) = spgu(ji,jj) + (  grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) )  &
-                           &                         + zld  * ( sshn     (ji+1,jj) - sshn     (ji,jj) )  ) * r1_e1u(ji,jj)
-                        spgv(ji,jj) = spgv(ji,jj) + (  grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) )  &
-                           &                         + zld  * ( sshn     (ji,jj+1) - sshn     (ji,jj) )  ) * r1_e2v(ji,jj)
+                        spgu(ji,jj) = spgu(ji,jj) + (  grav * ( pot_astro(ji+1,jj)     - pot_astro(ji,jj)     )  &
+                           &                         + zld  * ( ssh      (ji+1,jj,Nnn) - ssh      (ji,jj,Nnn) )  ) * r1_e1u(ji,jj)
+                        spgv(ji,jj) = spgv(ji,jj) + (  grav * ( pot_astro(ji,jj+1)     - pot_astro(ji,jj)     )  &
+                           &                         + zld  * ( ssh      (ji,jj+1,Nnn) - ssh      (ji,jj,Nnn) )  ) * r1_e2v(ji,jj)
                      END DO 
                   END DO

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynspg_exp.F90

@@ -49 +49 @@
       !!              momentum trend the surface pressure gradient :
       !!                      (ua,va) = (ua,va) + (spgu,spgv)
-      !!              where spgu = -1/rho0 d/dx(ps) = -g/e1u di( sshn )
-      !!                    spgv = -1/rho0 d/dy(ps) = -g/e2v dj( sshn )
+      !!              where spgu = -1/rho0 d/dx(ps) = -g/e1u di( ssh(Nnn) )
+      !!                    spgv = -1/rho0 d/dy(ps) = -g/e2v dj( ssh(Nnn) )
       !!
       !! ** Action :   (ua,va)   trend of horizontal velocity increased by 
@@ -74 +74 @@
          DO jj = 2, jpjm1                    ! now surface pressure gradient
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) * r1_e1u(ji,jj)
-               spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) * r1_e2v(ji,jj)
+               spgu(ji,jj) = - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)
+               spgv(ji,jj) = - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
             END DO 
          END DO

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynspg_ts.F90

@@ -130 +130 @@
       !!
       !! ** Action :
-      !!      -Update the filtered free surface at step "n+1"      : ssha
+      !!      -Update the filtered free surface at step "n+1"      : ssh(Naa)
       !!      -Update filtered barotropic velocities at step "n+1" : ua_b, va_b
       !!      -Compute barotropic advective fluxes at step "n"     : un_adv, vn_adv
@@ -440 +440 @@
             DO jj = 2, jpjm1
                DO ji = 2, jpim1 
-                  ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji+1,jj) ) >                &
-                     &      MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            &
-                     &      MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji+1,jj) + ht_0(ji+1,jj) )  &
-                     &                                                         > rn_wdmin1 + rn_wdmin2
-                  ll_tmp2 = ( ABS( sshn(ji+1,jj)            -  sshn(ji  ,jj))  > 1.E-12 ).AND.( &
-                     &      MAX(   sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                &
-                     &      MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
+                  ll_tmp1 = MIN(  ssh(ji,jj,Nnn)               ,  ssh(ji+1,jj,Nnn)                 )   >                       &
+                     &      MAX(                 - ht_0(ji,jj) ,                   - ht_0(ji+1,jj) ) .AND.                     &
+                     &      MAX(  ssh(ji,jj,Nnn) + ht_0(ji,jj) ,  ssh(ji+1,jj,Nnn) + ht_0(ji+1,jj) )   >   rn_wdmin1 + rn_wdmin2
+                     !
+                  ll_tmp2 = ABS(  ssh(ji+1,jj,Nnn)   -   ssh(ji,jj,Nnn)  )  > 1.E-12    .AND.    &
+                     &      MAX(  ssh(ji+1,jj,Nnn)   ,   ssh(ji,jj,Nnn)  )  >                    &
+                     &      MAX(-ht_0(ji+1,jj)       , -ht_0(ji,jj)      ) + rn_wdmin1 + rn_wdmin2
                   IF(ll_tmp1) THEN
                      zcpx(ji,jj) = 1.0_wp
                   ELSEIF(ll_tmp2) THEN
-                     ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
-                     zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &
-                                 &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) )
+                     ! no worries about  ssh(ji+1,jj,Nnn) -  ssh(ji,jj,Nnn) = 0, it won't happen ! here
+                     zcpx(ji,jj) = ABS(  ( ssh(ji+1,jj,Nnn) + ht_0(ji+1,jj) - ssh(ji,jj,Nnn) - ht_0(ji,jj) )   &
+                                 &     / ( ssh(ji+1,jj,Nnn)                 - ssh(ji,jj,Nnn)               )   )
                      zcpx(ji,jj) = MAX(  0._wp , MIN( zcpx(ji,jj) , 1._wp )  )
                   ELSE
@@ -458 +458 @@
                   ENDIF
                   !
-                  ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji,jj+1) ) >                &
-                     &      MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            &
-                     &      MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji,jj+1) + ht_0(ji,jj+1) )  &
-                     &                                                       > rn_wdmin1 + rn_wdmin2
-                  ll_tmp2 = ( ABS( sshn(ji,jj)              -  sshn(ji,jj+1))  > 1.E-12 ).AND.( &
-                     &      MAX(   sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                &
-                     &      MAX(  -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
+                  ll_tmp1 = MIN(  ssh(ji,jj,Nnn)               ,  ssh(ji,jj+1,Nnn)                 )   >                     &
+                     &      MAX(                 - ht_0(ji,jj) ,                   - ht_0(ji,jj+1) ) .AND.                   &
+                     &      MAX(  ssh(ji,jj,Nnn) + ht_0(ji,jj) ,  ssh(ji,jj+1,Nnn) + ht_0(ji,jj+1) )   > rn_wdmin1 + rn_wdmin2
+                     !
+                  ll_tmp2 = ABS(  ssh(ji,jj,Nnn)  -  ssh(ji,jj+1,Nnn) )  > 1.E-12   .AND.       &
+                     &      MAX(  ssh(ji,jj,Nnn)  ,  ssh(ji,jj+1,Nnn) )  >                      &
+                     &    ( MAX(-ht_0(ji,jj)      ,-ht_0(ji,jj+1)     ) + rn_wdmin1 + rn_wdmin2 )
   
                   IF(ll_tmp1) THEN
                      zcpy(ji,jj) = 1.0_wp
                   ELSE IF(ll_tmp2) THEN
-                     ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
-                     zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &
-                        &             / (sshn(ji,jj+1) - sshn(ji,jj  )) )
+                     ! no worries about   ssh(ji,jj+1,Nnn) -  ssh(ji,jj  ,Nnn) = 0, it won't happen ! here
+                     zcpy(ji,jj) = ABS( ( ssh(ji,jj+1,Nnn) + ht_0(ji,jj+1) - ssh(ji,jj,Nnn) - ht_0(ji,jj) )  &
+                        &             / ( ssh(ji,jj+1,Nnn)                 - ssh(ji,jj,Nnn)               )  )
                      zcpy(ji,jj) = MAX(  0._wp , MIN( zcpy(ji,jj) , 1.0_wp )  )
                   ELSE
@@ -481 +481 @@
             DO jj = 2, jpjm1
                DO ji = 2, jpim1
-                  zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   &
+                  zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) )   &
                      &                          * r1_e1u(ji,jj) * zcpx(ji,jj)  * wdrampu(ji,jj)  !jth
-                  zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   &
+                  zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) )   &
                      &                          * r1_e2v(ji,jj) * zcpy(ji,jj)  * wdrampv(ji,jj)  !jth
                END DO
@@ -492 +492 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) * r1_e1u(ji,jj)
-                  zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) * r1_e2v(ji,jj) 
+                  zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn)  ) * r1_e1u(ji,jj)
+                  zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn)  ) * r1_e2v(ji,jj) 
                END DO
             END DO
@@ -665 +665 @@
       !
       IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                    
-         sshn_e(:,:) =    sshn(:,:)            
+         sshn_e(:,:) =    ssh (:,:,Nnn)            
          un_e  (:,:) =    un_b(:,:)            
          vn_e  (:,:) =    vn_b(:,:)
@@ -674 +674 @@
          hvr_e (:,:) = r1_hv_n(:,:)
       ELSE                                ! CENTRED integration: start from BEFORE fields
-         sshn_e(:,:) =    sshb(:,:)
+         sshn_e(:,:) =    ssh (:,:,Nbb)
          un_e  (:,:) =    ub_b(:,:)         
          vn_e  (:,:) =    vb_b(:,:)
@@ -687 +687 @@
       !
       ! Initialize sums:
-      ua_b  (:,:) = 0._wp       ! After barotropic velocities (or transport if flux form)          
-      va_b  (:,:) = 0._wp
-      ssha  (:,:) = 0._wp       ! Sum for after averaged sea level
-      un_adv(:,:) = 0._wp       ! Sum for now transport issued from ts loop
-      vn_adv(:,:) = 0._wp
+      ua_b  (:,:)     = 0._wp       ! After barotropic velocities (or transport if flux form)          
+      va_b  (:,:)     = 0._wp
+      ssh   (:,:,Naa) = 0._wp       ! Sum for after averaged sea level
+      un_adv(:,:)     = 0._wp       ! Sum for now transport issued from ts loop
+      vn_adv(:,:)     = 0._wp
       !
       IF( ln_wd_dl ) THEN
@@ -1185 +1185 @@
          ENDIF
          !                                          ! Sum sea level
-         ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
+         ssh(:,:,Naa) = ssh(:,:,Naa) + za1 * ssha_e(:,:)
 
          !                                                 ! ==================== !
@@ -1223 +1223 @@
          END DO
       ELSE
-         ! At this stage, ssha has been corrected: compute new depths at velocity points
+         ! At this stage, ssh(Naa) has been corrected: compute new depths at velocity points
+!!gm KE conserving expression in Vector form 
+!         DO jj = 1, jpjm1
+!            DO ji = 1, jpim1      ! NO Vector Opt.
+!               zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) * (   e1e2t(ji  ,jj) * ssh(ji  ,jj,Naa)   &
+!                  &                                                         + e1e2t(ji+1,jj) * ssh(ji+1,jj,Naa)   )
+!               zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) * (   e1e2t(ji,jj  ) * ssh(ji,jj  ,Naa)   &
+!                  &                                                         + e1e2t(ji,jj+1) * ssh(ji,jj+1,Naa) )
+!            END DO
+!         END DO
+!! replace by the KE conserving expression in flux form
          DO jj = 1, jpjm1
             DO ji = 1, jpim1      ! NO Vector Opt.
-               zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj) * (   e1e2t(ji  ,jj) * ssha(ji  ,jj)   &
-                  &                                                          + e1e2t(ji+1,jj) * ssha(ji+1,jj)   )
-               zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj)  * r1_e1e2v(ji,jj) * (   e1e2t(ji,jj  ) * ssha(ji,jj  )   &
-                  &                                                          + e1e2t(ji,jj+1) * ssha(ji,jj+1) )
-            END DO
-         END DO
+               zsshu_a(ji,jj) = r1_2 * ( ssh(ji,jj,Naa) + ssh(ji+1,jj,Naa)  ) * ssumask(ji,jj)
+               zsshv_a(ji,jj) = r1_2 * ( ssh(ji,jj,Naa) + ssh(ji,jj+1,Naa)  ) * ssvmask(ji,jj)
+            END DO
+         END DO
+!!gm end
          CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
          !
-         DO jk=1,jpkm1
+         DO jk = 1, jpkm1
             ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * r1_Dt
             va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * r1_Dt

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/sshwzv.F90

@@ -56 +56 @@
       !!                ***  ROUTINE ssh_nxt  ***
       !!                   
-      !! ** Purpose :   compute the after ssh (ssha)
+      !! ** Purpose :   compute the after ssh(Naa)
       !!
       !! ** Method  : - Using the incompressibility hypothesis, the ssh increment
@@ -62 +62 @@
       !!      by the time step.
       !!
-      !! ** action  :   ssha, after sea surface height
+      !! ** action  :    ssh(Naa), after sea surface height
       !!
       !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
@@ -87 +87 @@
       !                                           !------------------------------!
 
-      IF(ln_wd_il)   CALL wad_lmt( sshb, z1_2rho0 * (emp_b(:,:) + emp(:,:)), rDt )
+      IF(ln_wd_il)   CALL wad_lmt( ssh(:,:,Nbb), z1_2rho0 * (emp_b(:,:) + emp(:,:)), rDt )
 
       CALL div_hor( kt )                               ! Horizontal divergence
@@ -99 +99 @@
       ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
       ! 
-      ssha(:,:) = (  sshb(:,:) - rDt * ( z1_2rho0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
+      ssh(:,:,Naa) = (  ssh(:,:,Nbb) - rDt * ( z1_2rho0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
       !
 #if defined key_agrif
@@ -107 +107 @@
       IF ( .NOT.ln_dynspg_ts ) THEN
          IF( ln_bdy ) THEN
-            CALL lbc_lnk( ssha, 'T', 1. )    ! Not sure that's necessary
-            CALL bdy_ssh( ssha )             ! Duplicate sea level across open boundaries
+            CALL lbc_lnk( ssh(:,:,Naa), 'T', 1. )    ! Not sure that's necessary
+            CALL bdy_ssh( ssh(:,:,Naa) )             ! Duplicate sea level across open boundaries
          ENDIF
       ENDIF
@@ -115 +115 @@
       !                                           !------------------------------!
       !
-      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask )
+      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssh(:,:,Naa), clinfo1=' ssha  - : ', mask1=tmask )
       !
       IF( ln_timing )   CALL timing_stop('ssh_nxt')
@@ -212 +212 @@
       !! ** Purpose :   achieve the sea surface  height time stepping by 
       !!              applying Asselin time filter and swapping the arrays
-      !!              ssha  already computed in ssh_nxt  
+      !!              ssh(Naa)  already computed in ssh_nxt  
       !!
       !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
-      !!              from the filter, see Leclair and Madec 2010) and swap :
-      !!                sshn = ssha + rn_atfp * ( sshb -2 sshn + ssha )
-      !!                            - rn_atfp * rn_Dt * ( emp_b - emp ) / rho0
-      !!                sshn = ssha
-      !!
-      !! ** action  : - sshb, sshn   : before & now sea surface height
-      !!                               ready for the next time step
+      !!              from the filter, see Leclair and Madec 2010) :
+      !!                ssh(Nnn) = ssh(Naa) + rn_atfp * ( ssh(Nbb) -2 ssh(Nnn) + ssh(Naa) )
+      !!                                  - rn_atfp * rn_Dt * ( emp_b - emp ) / rho0
+      !!              - swap the time-level indexes
+      !!
+      !! ** action  : - ssh(Nbb), ssh(Nnn) : new before & now sea surface height
+      !!                                     ready for the next time step
       !!
       !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
@@ -227 +227 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !
+      INTEGER  ::   isave   ! local integer
       REAL(wp) ::   zcoef   ! local scalar
       !!----------------------------------------------------------------------
@@ -240 +241 @@
       IF ( l_1st_euler ) THEN    !==  Euler time-stepping  ==!   no filter, just swap
          !
-         sshn(:,:) = ssha(:,:)               ! now    <-- after  (before already = now)
+         isave = Nnn
+         Nnn   = Naa                   ! now   <-- after  (before remains unchanged)
+         Naa   = isave                 ! after <-- previously now index
          !
       ELSE                       !==  Leap-Frog time-stepping  ==!   Asselin filter + swap
          !
          !                                   ! before <-- now filtered
-         sshb(:,:) = sshn(:,:) + rn_atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )
+         ssh(:,:,Nbb) = ssh(:,:,Nnn) + rn_atfp * ( ssh(:,:,Nbb) - 2 * ssh(:,:,Nnn) + ssh(:,:,Naa) )
          IF( .NOT.ln_linssh ) THEN           ! before <-- with forcing removed
             zcoef = rn_atfp * rn_Dt * r1_rho0
-            sshb(:,:) = sshb(:,:) - zcoef * (     emp_b(:,:) - emp   (:,:)   &
-               &                             -    rnf_b(:,:) + rnf   (:,:)   &
-               &                             + fwfisf_b(:,:) - fwfisf(:,:)   ) * ssmask(:,:)
+            ssh(:,:,Nbb) = ssh(:,:,Nbb) - zcoef * (     emp_b(:,:) - emp   (:,:)   &
+               &                                   -    rnf_b(:,:) + rnf   (:,:)   &
+               &                                   + fwfisf_b(:,:) - fwfisf(:,:)   ) * ssmask(:,:)
          ENDIF
-         sshn(:,:) = ssha(:,:)               ! now <-- after
-      ENDIF
-      !
-      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask )
+         isave = Nnn
+         Nnn   = Naa                   ! now    <-- after
+         Naa   = isave                 ! after  <-- previously now index
+      ENDIF
+      !
+      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssh(:,:,Nbb), clinfo1=' sshb  - : ', mask1=tmask )
       !
       IF( ln_timing )   CALL timing_stop('ssh_swp')

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/wet_dry.F90

@@ -117 +117 @@
 
 
-   SUBROUTINE wad_lmt( sshb1, sshemp, p2dt )
+   SUBROUTINE wad_lmt( pssh, sshemp, p2dt )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE wad_lmt  ***
@@ -127 +127 @@
       !! ** Action  : - calculate flux limiter and W/D flag
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   sshb1        !!gm DOCTOR names: should start with p !
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pssh      ! before sea-surface height
       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   sshemp
       REAL(wp)                , INTENT(in   ) ::   p2dt
@@ -178 +178 @@
                &         + MIN( zflxv(ji,jj) , 0._wp ) - MAX( zflxv(ji,  jj-1) , 0._wp ) 
             !
-            zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1
+            zdep2 = ht_0(ji,jj) + pssh(ji,jj) - rn_wdmin1
             IF( zdep2 <= 0._wp ) THEN     ! add more safty, but not necessary
-               sshb1(ji,jj) = rn_wdmin1 - ht_0(ji,jj)
+               pssh(ji,jj) = rn_wdmin1 - ht_0(ji,jj)
                IF(zflxu(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = 0._wp
                IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp
@@ -191 +191 @@
       !
       !           ! HPG limiter from jholt
-      wdramp(:,:) = min((ht_0(:,:) + sshb1(:,:) - rn_wdmin1)/(rn_wdmin0 - rn_wdmin1),1.0_wp)
+      wdramp(:,:) = min((ht_0(:,:) + pssh(:,:) - rn_wdmin1)/(rn_wdmin0 - rn_wdmin1),1.0_wp)
       !jth assume don't need a lbc_lnk here
       DO jj = 1, jpjm1
@@ -221 +221 @@
                !
                zdep1 = (zzflxp + zzflxn) * p2dt / ztmp
-               zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - p2dt * sshemp(ji,jj)
+               zdep2 = ht_0(ji,jj) + pssh(ji,jj) - rn_wdmin1 - p2dt * sshemp(ji,jj)
                !
                IF( zdep1 > zdep2 ) THEN

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/IOM/iom.F90

@@ -408 +408 @@
 
         i = 0
-        i = i + 1; fields(i)%vname="rdt";            fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="un";             fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="ub";             fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="vn";             fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="vb";             fields(i)%grid="grid_N_3D"  
-        i = i + 1; fields(i)%vname="tn";             fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="tb";             fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="sn";             fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="sb";             fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="sshn";           fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="sshb";           fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="rhop";           fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="kt";             fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="ndastp";         fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="adatrj";         fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="utau_b";         fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="vtau_b";         fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="qns_b";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="emp_b";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="sfx_b";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="en" ;            fields(i)%grid="grid_N_3D" 
-        i = i + 1; fields(i)%vname="avt_k";            fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="avm_k";            fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="dissl";          fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="sbc_hc_b";       fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="sbc_sc_b";       fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="qsr_hc_b";       fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="fraqsr_1lev";    fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="greenland_icesheet_mass"
-                                               fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="greenland_icesheet_timelapsed"
-                                               fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="greenland_icesheet_mass_roc"
-                                               fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="antarctica_icesheet_mass"
-                                               fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="antarctica_icesheet_timelapsed"
-                                               fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="antarctica_icesheet_mass_roc"
-                                               fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="frc_v";          fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="frc_t";          fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="frc_s";          fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="frc_wn_t";       fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="frc_wn_s";       fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="ssh_ini";        fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="e3t_ini";        fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="hc_loc_ini";     fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="sc_loc_ini";     fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="ssh_hc_loc_ini"; fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="ssh_sc_loc_ini"; fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="tilde_e3t_b";    fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="tilde_e3t_n";    fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="hdiv_lf";        fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="ub2_b";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="vb2_b";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="sshbb_e";        fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="ubb_e";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="vbb_e";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="sshb_e";         fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="ub_e";           fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="vb_e";           fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="fwf_isf_b";      fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="isf_sc_b";       fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="isf_hc_b";       fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="ssh_ibb";        fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="rnf_b";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="rnf_hc_b";       fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="rnf_sc_b";       fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="nn_fsbc";        fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="ssu_m";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="ssv_m";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="sst_m";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="sss_m";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="ssh_m";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="e3t_m";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="frq_m";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="avmb";           fields(i)%grid="grid_vector"
-        i = i + 1; fields(i)%vname="avtb";           fields(i)%grid="grid_vector"
-        i = i + 1; fields(i)%vname="ub2_i_b";        fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="vb2_i_b";        fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="ntime";          fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="Dsst";           fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="tmask";          fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="umask";          fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="vmask";          fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="smask";          fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="gdepw_n";        fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="e3t_n";          fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="e3u_n";          fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="e3v_n";          fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="surf_ini";       fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="e3t_b";          fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="hmxl_n";         fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="un_bf";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="vn_bf";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="hbl";            fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="hbli";           fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="wn";             fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="rdt"                             ;   fields(i)%grid="grid_scalar"
+        i = i + 1; fields(i)%vname="un"                              ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="ub"                              ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="vn"                              ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="vb"                              ;   fields(i)%grid="grid_N_3D"  
+        i = i + 1; fields(i)%vname="tn"                              ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="tb"                              ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="sn"                              ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="sb"                              ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname='sshn'                            ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname='sshb'                            ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="rhop"                            ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="kt"                              ;   fields(i)%grid="grid_scalar"
+        i = i + 1; fields(i)%vname="ndastp"                          ;   fields(i)%grid="grid_scalar"
+        i = i + 1; fields(i)%vname="adatrj"                          ;   fields(i)%grid="grid_scalar"
+        i = i + 1; fields(i)%vname="utau_b"                          ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="vtau_b"                          ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="qns_b"                           ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="emp_b"                           ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="sfx_b"                           ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="en"                              ;   fields(i)%grid="grid_N_3D" 
+        i = i + 1; fields(i)%vname="avt_k"                           ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="avm_k"                           ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="dissl"                           ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="sbc_hc_b"                        ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="sbc_sc_b"                        ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="qsr_hc_b"                        ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="fraqsr_1lev"                     ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="greenland_icesheet_mass"         ;   fields(i)%grid="grid_scalar"
+        i = i + 1; fields(i)%vname="greenland_icesheet_timelapsed"   ;   fields(i)%grid="grid_scalar"
+        i = i + 1; fields(i)%vname="greenland_icesheet_mass_roc"     ;   fields(i)%grid="grid_scalar"
+        i = i + 1; fields(i)%vname="antarctica_icesheet_mass"        ;   fields(i)%grid="grid_scalar"
+        i = i + 1; fields(i)%vname="antarctica_icesheet_timelapsed"  ;   fields(i)%grid="grid_scalar"
+        i = i + 1; fields(i)%vname="antarctica_icesheet_mass_roc"    ;   fields(i)%grid="grid_scalar"
+        i = i + 1; fields(i)%vname="frc_v"                           ;   fields(i)%grid="grid_scalar"
+        i = i + 1; fields(i)%vname="frc_t"                           ;   fields(i)%grid="grid_scalar"
+        i = i + 1; fields(i)%vname="frc_s"                           ;   fields(i)%grid="grid_scalar"
+        i = i + 1; fields(i)%vname="frc_wn_t"                        ;   fields(i)%grid="grid_scalar"
+        i = i + 1; fields(i)%vname="frc_wn_s"                        ;   fields(i)%grid="grid_scalar"
+        i = i + 1; fields(i)%vname="ssh_ini"                         ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="e3t_ini"                         ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="hc_loc_ini"                      ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="sc_loc_ini"                      ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="ssh_hc_loc_ini"                  ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="ssh_sc_loc_ini"                  ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="tilde_e3t_b"                     ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="tilde_e3t_n"                     ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="hdiv_lf"                         ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="ub2_b"                           ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="vb2_b"                           ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="sshbb_e"                         ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="ubb_e"                           ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="vbb_e"                           ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="sshb_e"                          ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="ub_e"                            ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="vb_e"                            ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="fwf_isf_b"                       ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="isf_sc_b"                        ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="isf_hc_b"                        ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="ssh_ibb"                         ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="rnf_b"                           ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="rnf_hc_b"                        ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="rnf_sc_b"                        ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="nn_fsbc"                         ;   fields(i)%grid="grid_scalar"
+        i = i + 1; fields(i)%vname="ssu_m"                           ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="ssv_m"                           ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="sst_m"                           ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="sss_m"                           ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="ssh_m"                           ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="e3t_m"                           ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="frq_m"                           ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="avmb"                            ;   fields(i)%grid="grid_vector"
+        i = i + 1; fields(i)%vname="avtb"                            ;   fields(i)%grid="grid_vector"
+        i = i + 1; fields(i)%vname="ub2_i_b"                         ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="vb2_i_b"                         ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="ntime"                           ;   fields(i)%grid="grid_scalar"
+        i = i + 1; fields(i)%vname="Dsst"                            ;   fields(i)%grid="grid_scalar"
+        i = i + 1; fields(i)%vname="tmask"                           ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="umask"                           ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="vmask"                           ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="smask"                           ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="gdepw_n"                         ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="e3t_n"                           ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="e3u_n"                           ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="e3v_n"                           ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="surf_ini"                        ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="e3t_b"                           ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="hmxl_n"                          ;   fields(i)%grid="grid_N_3D"
+        i = i + 1; fields(i)%vname="un_bf"                           ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="vn_bf"                           ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="hbl"                             ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="hbli"                            ;   fields(i)%grid="grid_N"
+        i = i + 1; fields(i)%vname="wn"                              ;   fields(i)%grid="grid_N_3D"
 
         IF( i-1 > max_rst_fields) THEN
@@ -523 +517 @@
       !!---------------------------------------------------------------------
 !sets enabled = .TRUE. for each field in restart file
-   CHARACTER(len=*) :: cdrst_file
+      CHARACTER(len=*) :: cdrst_file
 #if defined key_iomput
-   TYPE(xios_field) :: field_hdl
-   TYPE(xios_file) :: file_hdl
-   TYPE(xios_filegroup) :: filegroup_hdl
-   INTEGER :: i
-   CHARACTER(lc)  ::   clpath
+      TYPE(xios_field) :: field_hdl
+      TYPE(xios_file) :: file_hdl
+      TYPE(xios_filegroup) :: filegroup_hdl
+      INTEGER :: i
+      CHARACTER(lc)  ::   clpath
+      !!---------------------------------------------------------------------
 
 !set name of the restart file and enable available fields
@@ -567 +562 @@
 #endif
    END SUBROUTINE iom_set_rstw_active
+
 
    SUBROUTINE iom_set_rst_context( ) 

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/OBS/diaobs.F90

@@ -498 +498 @@
       USE dom_oce, ONLY : gdept_n, gdept_1d   ! Ocean space and time domain variables
       USE phycst , ONLY : rday                ! Physical constants
-      USE oce    , ONLY : tsn, un, vn, sshn   ! Ocean dynamics and tracers variables
+      USE oce    , ONLY : tsn, un, vn, ssh    ! Ocean dynamics and tracers variables
       USE phycst , ONLY : rday                ! Physical constants
 #if defined  key_si3
@@ -596 +596 @@
                zsurfvar(:,:) = tsn(:,:,1,jp_tem)
             CASE('sla')
-               zsurfvar(:,:) = sshn(:,:)
+               zsurfvar(:,:) = ssh(:,:,Nnn)
             CASE('sss')
                zsurfvar(:,:) = tsn(:,:,1,jp_sal)

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/OBS/obs_read_altbias.F90

@@ -28 +28 @@
       & e2t,   &
       & gphit
-   USE oce, ONLY : &           ! Model variables
-      & sshn
    USE obs_inter_h2d
    USE obs_utils               ! Various observation tools
@@ -35 +33 @@
 
    IMPLICIT NONE
-
-   !! * Routine accessibility
    PRIVATE
 
-   PUBLIC obs_rea_altbias     ! Read the altimeter bias
+   PUBLIC   obs_rea_altbias   ! Read the altimeter bias
 
    !!----------------------------------------------------------------------

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/OBS/obs_readmdt.F90

@@ -25 +25 @@
       &                    tmask, tmask_i, e1e2t, gphit, glamt
    USE obs_const, ONLY :   obfillflt      ! Fillvalue
-   USE oce      , ONLY :   sshn           ! Model variables
+   USE oce      , ONLY :   ssh            ! Model variables
 
    IMPLICIT NONE
@@ -216 +216 @@
           zarea = zarea + zdxdy
           zeta1 = zeta1 + mdt(ji,jj) * zdxdy
-          zeta2 = zeta2 + sshn (ji,jj) * zdxdy
+          zeta2 = zeta2 + ssh(ji,jj,Nnn) * zdxdy
         END DO      
       END DO
@@ -241 +241 @@
       IF(lwp) THEN
          WRITE(numout,*)
-         WRITE(numout,*) ' obs_readmdt : rn_mdtcutoff     = ', rn_mdtcutoff
+         WRITE(numout,*) ' obs_readmdt : rn_mdtcutoff  = ', rn_mdtcutoff
          WRITE(numout,*) ' -----------   zcorr_mdt     = ', zcorr_mdt
          WRITE(numout,*) '               zcorr_bcketa  = ', zcorr_bcketa
          WRITE(numout,*) '               zcorr         = ', zcorr
-         WRITE(numout,*) '               nn_msshc        = ', nn_msshc
+         WRITE(numout,*) '               nn_msshc      = ', nn_msshc
       ENDIF
 

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/OBS/obs_sstbias.F90

@@ -27 +27 @@
       & gphit, &
       & glamt
-   USE oce, ONLY : &           ! Model variables
-      & sshn
    USE obs_inter_h2d
    USE obs_utils               ! Various observation tools

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/SBC/sbccpl.F90

@@ -32 +32 @@
    USE cpl_oasis3     ! OASIS3 coupling
    USE geo2ocean      ! 
-   USE oce     , ONLY : tsn, un, vn, sshn, ub, vb, sshb, fraqsr_1lev
+   USE oce     , ONLY : tsn, un, vn, ssh, ub, vb, fraqsr_1lev
    USE ocealb         ! 
    USE eosbn2         ! 
@@ -2467 +2467 @@
          IF( ln_apr_dyn ) THEN  
             IF( kt /= nit000 ) THEN  
-               ztmp1(:,:) = sshb(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )  
+               ztmp1(:,:) = ssh(:,:,Nbb) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )  
             ELSE  
-               ztmp1(:,:) = sshb(:,:)  
+               ztmp1(:,:) = ssh(:,:,Nbb)  
             ENDIF  
          ELSE  
-            ztmp1(:,:) = sshn(:,:)  
+            ztmp1(:,:) = ssh(:,:,Nnn)  
          ENDIF  
          CALL cpl_snd( jps_wlev  , isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info ) 
@@ -2482 +2482 @@
          !                          ! removed inverse barometer ssh when Patm
          !                          forcing is used (for sea-ice dynamics)
-         IF( ln_apr_dyn ) THEN   ;   ztmp1(:,:) = sshb(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
-         ELSE                    ;   ztmp1(:,:) = sshn(:,:)
+         IF( ln_apr_dyn ) THEN   ;   ztmp1(:,:) = ssh(:,:,Nbb) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
+         ELSE                    ;   ztmp1(:,:) = ssh(:,:,Nnn)
          ENDIF
          CALL cpl_snd( jps_ssh   , isec, RESHAPE ( ztmp1            , (/jpi,jpj,1/) ), info )
-
       ENDIF
       !                                                        ! SSS

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/SBC/sbcfwb.F90

@@ -127 +127 @@
             a_fwb_b = a_fwb                           ! mean sea level taking into account the ice+snow
                                                       ! sum over the global domain
-            a_fwb   = glob_sum( e1e2t(:,:) * ( sshn(:,:) + snwice_mass(:,:) * r1_rho0 ) )
+            a_fwb   = glob_sum( e1e2t(:,:) * ( ssh(:,:,Nnn) + snwice_mass(:,:) * r1_rho0 ) )
             a_fwb   = a_fwb * 1.e+3 / ( area * rday * 365. )     ! convert in Kg/m3/s = mm/s
 !!gm        !                                                      !!bug 365d year 

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/SBC/sbcice_cice.F90

@@ -227 +227 @@
       IF( .NOT.ln_rstart ) THEN
          IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area
-            sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rho0
-            sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rho0
-
-!!gm This should be put elsewhere....   (same remark for limsbc)
-!!gm especially here it is assumed zstar coordinate, but it can be ztilde....
-            IF( .NOT.ln_linssh ) THEN
+            ssh(:,:,Nnn) = ssh(:,:,Nnn) - snwice_mass(:,:) * r1_rho0
+            ssh(:,:,Nbb) = ssh(:,:,Nbb) - snwice_mass(:,:) * r1_rho0
+            !
+            IF( .NOT.ln_linssh ) THEN     ! modified the now and before vertical mesh and scale factors 
                !
-               DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
-                  e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) * r1_ht_0(:,:) )
-                  e3t_b(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshb(:,:) * r1_ht_0(:,:) )
-               END DO
-               e3t_a(:,:,:) = e3t_b(:,:,:)
-               ! Reconstruction of all vertical scale factors at now and before time-steps
-               ! =============================================================================
-               ! Horizontal scale factor interpolations
-               ! --------------------------------------
-               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
-               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
-               CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
-               CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
-               CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
-               ! Vertical scale factor interpolations
-               ! ------------------------------------
-               CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )
-               CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
-               CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
-               CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
-               CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
-               ! t- and w- points depth
-               ! ----------------------
-               gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
-               gdepw_n(:,:,1) = 0.0_wp
-               gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
-               DO jk = 2, jpk
-                  gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk)
-                  gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1)
-                  gde3w_n(:,:,jk) = gdept_n(:,:,jk  ) - sshn   (:,:)
-               END DO
+               !                             !* BEFORE fields : 
+               CALL ssh2e3_before               ! set:      hu , hv , r1_hu, r1_hv 
+               !                                    !  e3t, e3w, e3u, e3uw, e3v, e3vw
+               !
+               !                             !* NOW fields : 
+               CALL ssh2e3_now                  ! set: ht , hu , hv , r1_hu, r1_hv
+               !                                !      e3t, e3w, e3u, e3uw, e3v, e3vw, e3f
+               !                                !      gdept_n, gdepw_n, gde3w_n
             ENDIF
          ENDIF
@@ -290 +265 @@
       ENDIF
 
-      ztmp(:,:)=0.0
+      ztmp(:,:) = 0._wp
 
 ! Aggregate ice concentration already set in cice_sbc_out (or cice_sbc_init on 

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/SBC/sbcssm.F90

@@ -75 +75 @@
          sss_m(:,:) = zts(:,:,jp_sal)
          !                          ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
-         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
-         ELSE                    ;   ssh_m(:,:) = sshn(:,:)
+         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = ssh(:,:,Nnn) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
+         ELSE                    ;   ssh_m(:,:) = ssh(:,:,Nnn)
          ENDIF
          !
@@ -98 +98 @@
             sss_m(:,:) = zcoef * zts(:,:,jp_sal)
             !                          ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
-            IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = zcoef * ( sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) )
-            ELSE                    ;   ssh_m(:,:) = zcoef * sshn(:,:)
+            IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = zcoef * ( ssh(:,:,Nnn) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) )
+            ELSE                    ;   ssh_m(:,:) = zcoef *   ssh(:,:,Nnn)
             ENDIF
             !
@@ -126 +126 @@
          sss_m(:,:) = sss_m(:,:) + zts(:,:,jp_sal)
          !                          ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
-         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
-         ELSE                    ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:)
+         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = ssh_m(:,:) + ssh(:,:,Nnn) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
+         ELSE                    ;   ssh_m(:,:) = ssh_m(:,:) + ssh(:,:,Nnn)
          ENDIF
          !
@@ -247 +247 @@
          ENDIF
          sss_m(:,:) = tsn  (:,:,1,jp_sal)
-         ssh_m(:,:) = sshn (:,:)
+         ssh_m(:,:) = ssh  (:,:,Nnn)
          e3t_m(:,:) = e3t_n(:,:,1)
          frq_m(:,:) = 1._wp

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/trasbc.F90

@@ -125 +125 @@
          DO ji = fs_2, fs_jpim1   ! vector opt.
             IF ( ll_wd ) THEN     ! If near WAD point limit the flux for now
-               IF ( sshn(ji,jj) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN
+               IF ( ssh(ji,jj,Nnn) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN
                   sbc_tsc(ji,jj,jp_tem) = r1_rho0_rcp * qns(ji,jj)   ! non solar heat flux
-               ELSE IF ( sshn(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN
+               ELSE IF ( ssh(ji,jj,Nnn) + ht_0(ji,jj) >  rn_wdmin1 ) THEN
                   sbc_tsc(ji,jj,jp_tem) = r1_rho0_rcp * qns(ji,jj) &
-                       &                * tanh ( 5._wp * ( ( sshn(ji,jj) + ht_0(ji,jj) -  rn_wdmin1 ) * r_rn_wdmin1 ) )
+                     &                  * tanh ( 5._wp * ( ( ssh(ji,jj,Nnn) + ht_0(ji,jj) - rn_wdmin1 ) * r_rn_wdmin1 ) )
                ELSE
                   sbc_tsc(ji,jj,jp_tem) = 0._wp

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRD/trdpen.F90

@@ -114 +114 @@
                                 !IF( ln_linssh ) THEN                   ! cst volume : ssh term (otherwise include in e3t variation)
                                 !   ALLOCATE( z2d(jpi,jpj) )
-                                !   z2d(:,:) = ( ssha(:,:) - sshb(:,:) )                 &
+                                !   z2d(:,:) = ( ssh(:,:,Naa) - ssh(:,:,Nbb) )          &
                                 !      &     * (   dPE_dt(:,:,1) * tsn(:,:,1,jp_tem)    &
                                 !      &         + dPE_ds(:,:,1) * tsn(:,:,1,jp_sal)  ) / ( e3t_n(:,:,1) * pdt )

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/nemogcm.F90

@@ -424 +424 @@
          !
       ENDIF
-      !                                      ! Domain decomposition
-      CALL mpp_init                          ! MPP
-
-      ! Now we know the dimensions of the grid and numout has been set: we can allocate arrays
-      CALL nemo_alloc()
+      
+      CALL mpp_init                          ! MPP domain decomposition
+
+      CALL dom_nam                           ! set the time (Nt) by reading namelist namrun and namdom
+
+      !                                      ! allocate arrays (here, as we know the dimensions of both 
+      !                                      !                  the grid and number of time-level used, and
+      CALL nemo_alloc()                      !                  numout has been set: we can allocate arrays)
 
       !                             !-------------------------------!

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/oce.F90

@@ -34 +34 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub_b   ,  un_b  ,  ua_b  !: Barotropic velocities at u-point [m/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vb_b   ,  vn_b  ,  va_b  !: Barotropic velocities at v-point [m/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshb   ,  sshn  , ssha   !: sea surface height at t-point [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ssh                    !: sea surface height at t-point    [m]
 
 !!gm??   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ssh  , sshu   , sshv   !: sea surface height at t-, u- v-points [m]
@@ -93 +93 @@
          &      rhd  (jpi,jpj,jpk)      , rhop (jpi,jpj,jpk)                              , STAT=ierr(1) )
          !
-      ALLOCATE( sshb  (jpi,jpj)  , sshn   (jpi,jpj) , ssha(jpi,jpj)   ,     &
-         &      ub_b(jpi,jpj)     , un_b(jpi,jpj)   , ua_b(jpi,jpj)   ,     &
+      ALLOCATE( ub_b(jpi,jpj)     , un_b(jpi,jpj)   , ua_b(jpi,jpj)   ,     &
          &      vb_b(jpi,jpj)     , vn_b(jpi,jpj)   , va_b(jpi,jpj)   ,     &
          &      spgu  (jpi,jpj)   , spgv(jpi,jpj)                     ,     &
@@ -102 +101 @@
          &      grui(jpi,jpj)     , grvi(jpi,jpj)                     ,     &
          &      riceload(jpi,jpj)                                     , STAT=ierr(2) )
+         !
+         
+         
+      ALLOCATE( ssh(jpi,jpj,Nt) , STAT=ierr(5) )
+
+
          !
       ALLOCATE( fraqsr_1lev(jpi,jpj) , STAT=ierr(3) )

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/par_oce.F90

@@ -13 +13 @@
    PUBLIC
 
-
-   !!
-   INTEGER, PUBLIC ::   Nnn, Np1, Nm1   ! =now, before, after
-
-   INTEGER, PUBLIC ::   Nb, Nn, Na      ! before, now, after  index
-
-
-
-   !!----------------------------------------------------------------------
-   !!                   namcfg namelist parameters
-   !!----------------------------------------------------------------------
-   LOGICAL       ::   ln_read_cfg      !: (=T) read the domain configuration file or (=F) not
-   CHARACTER(lc) ::      cn_domcfg        !: filename the configuration file to be read
-   LOGICAL       ::   ln_write_cfg     !: (=T) create the domain configuration file
-   CHARACTER(lc) ::      cn_domcfg_out    !: filename the configuration file to be read
-   !
-   LOGICAL       ::   ln_use_jattr     !: input file read offset
-   !                                   !  Use file global attribute: open_ocean_jstart to determine start j-row 
-   !                                   !  when reading input from those netcdf files that have the 
-   !                                   !  attribute defined. This is designed to enable input files associated 
-   !                                   !  with the extended grids used in the under ice shelf configurations to 
-   !                                   !  be used without redundant rows when the ice shelves are not in use.
-   ! 
-
    !!---------------------------------------------------------------------
    !! Domain Matrix size 
    !!---------------------------------------------------------------------
-   ! configuration name & resolution   (required only in ORCA family case)
-   CHARACTER(lc) ::   cn_cfg           !: name of the configuration
-   INTEGER       ::   nn_cfg           !: resolution of the configuration 
+   ! time dimension and index
+   INTEGER, PUBLIC ::   Nt             !: number of time-level used
+   INTEGER, PUBLIC ::   Nbb, Nnn, Naa  !: before, now, after time-level indices
 
    ! global domain size               !!! * total computational domain *
@@ -48 +24 @@
    INTEGER       ::   jpjglo           !: 2nd    -                  -    --> j-direction
    INTEGER       ::   jpkglo           !: 3nd    -                  -    --> k levels
+!!gm  to be used in futur (?):
+!!   INTEGER, PUBLIC ::   Niglo , Njglo , Nkglo !: global domain size
 
    ! global domain size for AGRIF     !!! * total AGRIF computational domain *
@@ -65 +43 @@
    INTEGER, PUBLIC ::   jpimax! = ( jpiglo-2*nn_hls + (jpni-1) ) / jpni + 2*nn_hls !: maximum jpi
    INTEGER, PUBLIC ::   jpjmax! = ( jpjglo-2*nn_hls + (jpnj-1) ) / jpnj + 2*nn_hls !: maximum jpj
+!!gm  to be used in futur (?):
+!!   INTEGER, PUBLIC ::   Ni  , Nj   , Nk     !: local domain size
 
    !!---------------------------------------------------------------------
@@ -85 +65 @@
 
    !!----------------------------------------------------------------------
+   !!                   namcfg namelist parameters
+   !!----------------------------------------------------------------------
+   LOGICAL       ::   ln_read_cfg      !: (=T) read the domain configuration file or (=F) not
+   CHARACTER(lc) ::      cn_domcfg        !: filename the configuration file to be read
+   LOGICAL       ::   ln_write_cfg     !: (=T) create the domain configuration file
+   CHARACTER(lc) ::      cn_domcfg_out    !: filename the configuration file to be read
+   LOGICAL       ::   ln_use_jattr     !: input file read offset
+   !                                   !  Use file global attribute: open_ocean_jstart to determine start j-row 
+   !                                   !  when reading input from those netcdf files that have the 
+   !                                   !  attribute defined. This is designed to enable input files associated 
+   !                                   !  with the extended grids used in the under ice shelf configurations to 
+   !                                   !  be used without redundant rows when the ice shelves are not in use.
+
+   ! configuration name & resolution   (required only in ORCA family case)
+   CHARACTER(lc) ::   cn_cfg           !: name of the configuration
+   INTEGER       ::   nn_cfg           !: resolution of the configuration 
+
+   !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
    !! $Id$ 

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/stpctl.F90

@@ -102 +102 @@
       !                                   !==  test of extrema  ==!
       IF( ll_wd ) THEN
-         zmax(1) = MAXVAL(  ABS( sshn(:,:) + ssh_ref*tmask(:,:,1) )  )        ! ssh max 
+         zmax(1) = MAXVAL(  ABS( ssh(:,:,Nnn) + ssh_ref*tmask(:,:,1) )  )        ! ssh max 
       ELSE
-         zmax(1) = MAXVAL(  ABS( sshn(:,:) )  )                               ! ssh max
+         zmax(1) = MAXVAL(  ABS( ssh(:,:,Nnn) )  )                               ! ssh max
       ENDIF
       zmax(2) = MAXVAL(  ABS( un(:,:,:) )  )                                  ! velocity max (zonal only)
@@ -129 +129 @@
          &  ISNAN( zmax(1) + zmax(2) + zmax(3) )  ) THEN   ! NaN encounter in the tests
          IF( lk_mpp ) THEN
-            CALL mpp_maxloc( ABS(sshn)        , ssmask(:,:)  , zzz, iih , ijh        )
+            CALL mpp_maxloc( ABS(ssh(:,:,Nnn)), ssmask(:,:)  , zzz, iih , ijh        )
             CALL mpp_maxloc( ABS(un)          , umask (:,:,:), zzz, iiu , iju , iku  )
             CALL mpp_minloc( tsn(:,:,:,jp_sal), tmask (:,:,:), zzz, iis1, ijs1, iks1 )
             CALL mpp_maxloc( tsn(:,:,:,jp_sal), tmask (:,:,:), zzz, iis2, ijs2, iks2 )
          ELSE
-            iloch  = MINLOC( ABS( sshn(:,:)   )                               )
-            ilocu  = MAXLOC( ABS( un  (:,:,:) )                               )
+            iloch  = MINLOC( ABS( ssh(:,:,Nnn) )                              )
+            ilocu  = MAXLOC( ABS( un (:,:,:)   )                              )
             ilocs1 = MINLOC( tsn(:,:,:,jp_sal) , mask = tmask(:,:,:) == 1._wp )
             ilocs2 = MAXLOC( tsn(:,:,:,jp_sal) , mask = tmask(:,:,:) == 1._wp )

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OFF/dtadyn.F90

@@ -142 +142 @@
          emp_b  (:,:)   = sf_dyn(jf_empb)%fnow(:,:,1) * tmask(:,:,1)    ! E-P
          zemp   (:,:)   = ( 0.5_wp * ( emp(:,:) + emp_b(:,:) ) + rnf(:,:) + fwbcorr ) * tmask(:,:,1)
-         CALL dta_dyn_ssh( kt, zhdivtr, sshb, zemp, ssha, e3t_a(:,:,:) )  !=  ssh, vertical scale factor & vertical transport
+         CALL dta_dyn_ssh( kt, zhdivtr, ssh(:,:,Nbb), zemp, ssh(:,:,Naa) )  !=  ssh, vertical scale factor & vertical transport
+!!
+!!gm BUG ?  ssh after computed but no swap so, not used in the restart....
+!!
          DEALLOCATE( zemp , zhdivtr )
          !                                           Write in the tracer restart file
@@ -148 +151 @@
          IF( lrst_trc ) THEN
             IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) 'dta_dyn_ssh : ssh field written in tracer restart file at it= ', kt,' date= ', ndastp
-            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-            CALL iom_rstput( kt, nitrst, numrtw, 'sshn', ssha )
-            CALL iom_rstput( kt, nitrst, numrtw, 'sshb', sshn )
+            IF(lwp) WRITE(numout,*) 'dta_dyn : ssh field written in tracer restart file at it= ', kt,' date= ', ndastp
+            IF(lwp) WRITE(numout,*) '~~~~~~~'
+            CALL iom_rstput( kt, nitrst, numrtw, 'sshn', ssh(:,:,Nnn) )
+            CALL iom_rstput( kt, nitrst, numrtw, 'sshb', ssh(:,:,Nbb) )
          ENDIF
       ENDIF
@@ -313 +316 @@
       !
       IF( .NOT.ln_linssh ) THEN
-        IF( .NOT. sf_dyn(jf_uwd)%ln_clim .AND. ln_rsttr .AND.    &                     ! Restart: read in restart file
-           iom_varid( numrtr, 'sshn', ldstop = .FALSE. ) > 0 ) THEN
-           IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation'
-           CALL iom_get( numrtr, jpdom_autoglo, 'sshn', sshn(:,:)   )
-           CALL iom_get( numrtr, jpdom_autoglo, 'sshb', sshb(:,:)   )
-        ELSE
-           IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation'
-           CALL iom_open( 'restart', inum )
-           CALL iom_get( inum, jpdom_autoglo, 'sshn', sshn(:,:)   )
-           CALL iom_get( inum, jpdom_autoglo, 'sshb', sshb(:,:)   )
-           CALL iom_close( inum )                                        ! close file
-        ENDIF
-        !
-        DO jk = 1, jpkm1
-           e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) )
-        ENDDO
-        e3t_a(:,:,jpk) = e3t_0(:,:,jpk)
-
-        ! Horizontal scale factor interpolations
-        ! --------------------------------------
-        CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
-        CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
-
-        ! Vertical scale factor interpolations
-        ! ------------------------------------
-        CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n(:,:,:), 'W' )
-  
-        e3t_b(:,:,:)  = e3t_n(:,:,:)
-        e3u_b(:,:,:)  = e3u_n(:,:,:)
-        e3v_b(:,:,:)  = e3v_n(:,:,:)
-
-        ! t- and w- points depth
-        ! ----------------------
-        gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
-        gdepw_n(:,:,1) = 0.0_wp
-
-        DO jk = 2, jpk
-           DO jj = 1,jpj
-              DO ji = 1,jpi
-                !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere
-                !    tmask = wmask, ie everywhere expect at jk = mikt
-                                                                   ! 1 for jk =
-                                                                   ! mikt
-                 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
-                 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
-                 gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
-                     &                + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk))
-              END DO
-           END DO
-        END DO
-
-        gdept_b(:,:,:) = gdept_n(:,:,:)
-        gdepw_b(:,:,:) = gdepw_n(:,:,:)
-        !
+         IF( .NOT. sf_dyn(jf_uwd)%ln_clim .AND. ln_rsttr .AND.    &                     ! Restart: read in restart file
+            iom_varid( numrtr, 'sshn', ldstop = .FALSE. ) > 0 ) THEN
+            IF(lwp) WRITE(numout,*) ' ssh forcing fields read in the restart file for initialisation'
+            CALL iom_get( numrtr, jpdom_autoglo, 'sshn', ssh(:,:,Nnn)   )
+            CALL iom_get( numrtr, jpdom_autoglo, 'sshb', ssh(:,:,Nbb)   )
+         ELSE
+            IF(lwp) WRITE(numout,*) ' ssh forcing fields read in the restart file for initialisation'
+            CALL iom_open( 'restart', inum )
+            CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh(:,:,Nnn)   )
+            CALL iom_get( inum, jpdom_autoglo, 'sshb', ssh(:,:,Nbb)   )
+            CALL iom_close( inum )                                        ! close file
+         ENDIF
+         !
+         !                    !== Set of all other vertical mesh fields  ==!  (now and before)      
+         !
+         !                          !* BEFORE fields : 
+         CALL ssh2e3_before               ! set:      hu , hv , r1_hu, r1_hv 
+         !                                !      e3t, e3w, e3u, e3uw, e3v, e3vw        (from 1 to jpkm1)
+         !
+         !                                ! set jpk level one to the e3._0 values
+         e3t_b(:,:,jpk) = e3t_0(:,:,jpk)  ;   e3u_b(:,:,jpk) =  e3w_0(:,:,jpk)  ;   e3v_b(:,:,jpk) =  e3v_0(:,:,jpk)
+         e3w_b(:,:,jpk) = e3w_0(:,:,jpk)  ;  e3uw_b(:,:,jpk) = e3uw_0(:,:,jpk)  ;  e3vw_b(:,:,jpk) = e3vw_0(:,:,jpk)
+         !
+         !                          !* NOW fields : 
+         CALL ssh2e3_now                  ! set: ht , hu , hv , r1_hu, r1_hv
+         !                                !      e3t, e3w, e3u, e3uw, e3v, e3vw, e3f   (from 1 to jpkm1)
+         !                                !      gdept_n, gdepw_n, gde3w_n
+!!gm issue?   gdept_n, gdepw_n, gde3w_n never defined at jpk
+         !
+         !                                ! set one for all last level to the e3._0 value
+         e3t_n(:,:,jpk) = e3t_0(:,:,jpk)  ;   e3u_n(:,:,jpk) =  e3w_0(:,:,jpk)  ;   e3v_n(:,:,jpk) =  e3v_0(:,:,jpk)
+         e3w_n(:,:,jpk) = e3w_0(:,:,jpk)  ;  e3uw_n(:,:,jpk) = e3uw_0(:,:,jpk)  ;  e3vw_n(:,:,jpk) = e3vw_0(:,:,jpk)
+         e3f_n(:,:,jpk) = e3f_0(:,:,jpk)
+         !
+         !                          !* AFTER fields : (last level for OPA, 3D required for AGRIF initialisation)
+         e3t_a(:,:,:) = e3t_n(:,:,:)   ;   e3u_a(:,:,:) = e3u_n(:,:,:)   ;   e3v_a(:,:,:) = e3v_n(:,:,:)
+         !
       ENDIF
       !
@@ -430 +416 @@
       INTEGER             :: ji, jj, jk
       REAL(wp)            :: zcoef
+      REAL(wp), DIMENSION(jpi,jpj) ::   zssht_h, zsshu_h, zsshv_h
       !!---------------------------------------------------------------------
 
@@ -438 +425 @@
       ENDIF
 
-      sshb(:,:) = sshn(:,:) + rn_atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )   ! before <-- now filtered
-      sshn(:,:) = ssha(:,:)
-
-      e3t_n(:,:,:) = e3t_a(:,:,:)
+      ssh(:,:,Nbb) = ssh(:,:,Nnn) + rn_atfp * ( ssh(:,:,Nbb) - 2 * ssh(:,:,Nnn) + ssh(:,:,Naa) )   ! before <-- now filtered
+      ssh(:,:,Nnn) = ssh(:,:,Naa)
 
       ! Reconstruction of all vertical scale factors at now and before time steps
       ! =============================================================================
-
-      ! Horizontal scale factor interpolations
-      ! --------------------------------------
-      CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
-      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
-
-      ! Vertical scale factor interpolations
-      ! ------------------------------------
-      CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' )
-
+      !
+      !                                   !==  now ssh  ==!  (u- and v-points)
+      DO jj = 2, jpjm1   ;   DO ji = 2, jpim1
+         zsshu_h(ji,jj) = 0.5_wp * ( ssh(ji,jj,Nnn) + ssh(ji+1,jj,Nnn) ) * ssumask(ji,jj)
+         zsshv_h(ji,jj) = 0.5_wp * ( ssh(ji,jj,Nnn) + ssh(ji,jj+1,Nnn) ) * ssvmask(ji,jj)
+      END DO             ;   END DO      
+      CALL lbc_lnk_multi( zsshu_h(:,:), 'U', 1._wp , zsshv_h(:,:), 'V', 1._wp )
+      !
+      !                                   !==  after depths and its inverse  ==! 
+         hu_n(:,:) = hu_0(:,:) + zsshu_h(:,:)
+         hv_n(:,:) = hv_0(:,:) + zsshv_h(:,:)
+      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )
+      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
+      !
+      !                                   !==  now scale factors  ==!  (e3t , e3u , e3v)
+      zssht_h(:,:) = ssh    (:,:,Nnn) * r1_ht_0(:,:)     ! t-point
+      zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)     ! u-point
+      zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)     ! v-point
+      DO jk = 1, jpkm1
+         e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
+         e3u_n(:,:,jk) = e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) )
+         e3v_n(:,:,jk) = e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) )
+         e3w_n(:,:,jk) = e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk) , tmask(:,:,jk+1) )
+      END DO
+      !
       e3t_b(:,:,:)  = e3t_n(:,:,:)
       e3u_b(:,:,:)  = e3u_n(:,:,:)
@@ -475 +475 @@
       END DO
       !
+      zssht_h(:,:) = 1._wp + zssht_h(:,:)               ! t-point
+      !
+      IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness 
+         DO jk = 1, jpkm1
+            gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:)
+            gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:)
+            gde3w_n(:,:,jk) =   gdept_n(:,:,jk) - ssh    (:,:,Nnn)
+         END DO
+      ELSE                    ! no ISF cavities 
+         DO jk = 1, jpkm1
+            gdept_n(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:)
+            gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:)
+            gde3w_n(:,:,jk) = gdept_n(:,:,jk) - ssh    (:,:,Nnn)
+         END DO
+      ENDIF
+      !
       gdept_b(:,:,:) = gdept_n(:,:,:)
       gdepw_b(:,:,:) = gdepw_n(:,:,:)
@@ -481 +497 @@
    
 
-   SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb,  pemp, pssha, pe3ta )
+   SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb,  pemp, pssha )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE dta_dyn_wzv  ***
@@ -502 +518 @@
       !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
       !!----------------------------------------------------------------------
-      INTEGER,                                   INTENT(in )    :: kt        !  time-step
-      REAL(wp), DIMENSION(jpi,jpj,jpk)          , INTENT(in )   :: phdivtr   ! horizontal divergence transport
-      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(in )   :: psshb     ! now ssh
-      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(in )   :: pemp      ! evaporation minus precipitation
+      INTEGER,                                    INTENT(in   )    :: kt        !  time-step
+      REAL(wp), DIMENSION(jpi,jpj,jpk)          , INTENT(in   )   :: phdivtr   ! horizontal divergence transport
+      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(in   )   :: psshb     ! now ssh
+      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(in   )   :: pemp      ! evaporation minus precipitation
       REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(inout) :: pssha     ! after ssh
-      REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL, INTENT(out)   :: pe3ta     ! after vertical scale factor
       !
       INTEGER                       :: jk
@@ -518 +533 @@
       END DO
       !                                                ! Sea surface  elevation time-stepping
-      pssha(:,:) = ( psshb(:,:) - rDt * ( r1_rho0 * pemp(:,:)  + zhdiv(:,:) ) ) * ssmask(:,:)
-      !                                                 ! 
-      !                                                 ! After acale factors at t-points ( z_star coordinate )
-      DO jk = 1, jpkm1
-        pe3ta(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + pssha(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) )
-      END DO
+      pssha(:,:) = ( psshb(:,:) - rDt * ( r1_rho0 * pemp(:,:) + zhdiv(:,:) ) ) * ssmask(:,:)
       !
    END SUBROUTINE dta_dyn_ssh

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OFF/nemogcm.F90

@@ -265 +265 @@
          !
       ENDIF
-      !                                      ! Domain decomposition
-      CALL mpp_init                          ! MPP
-
-      ! Now we know the dimensions of the grid and numout has been set: we can allocate arrays
-      CALL nemo_alloc()
+      CALL mpp_init                          ! MPP domain decomposition
+
+      CALL dom_nam                           ! set the time (Nt) by reading namelist namrun and namdom
+
+      !                                      ! allocate arrays (here, as we know the dimensions of both 
+      !                                      !                  the grid and number of time-level used, and
+      CALL nemo_alloc()                      !                  numout has been set: we can allocate arrays)
 
       !                             !-------------------------------!

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/SAO/sao_read.F90

@@ -9 +9 @@
    USE par_kind, ONLY: lc
    USE netcdf
-   USE oce,     ONLY: tsn, sshn
+   USE oce,     ONLY: tsn, ssh
    USE dom_oce, ONLY: nlci, nlcj, nimpp, njmpp, tmask
    USE par_oce, ONLY: jpi, jpj, jpk
@@ -77 +77 @@
 
       IF (TRIM(filename) == 'nofile') THEN
-         tsn (:,:,:,:) = fbrmdi
-         sshn(:,:)     = fbrmdi
+         tsn(:,:,:,:) = fbrmdi
+         ssh(:,:,Nnn) = fbrmdi
       ELSE
          WRITE(numout,*) "Opening :", TRIM(filename)
@@ -133 +133 @@
             WHERE(temp_sshn(:,:) == fill_val) temp_sshn(:,:) = fbrmdi
 
-            ! Initialise tsn, sshn to fbrmdi
+            ! Initialise tsn, ssh to fbrmdi
             tsn(:,:,:,:) = fbrmdi
-            sshn(:,:) = fbrmdi
+            ssh(:,:,Nnn) = fbrmdi
 
             ! Mask out missing data index
             tsn(1:nlci,1:nlcj,1:jpk,1) = temp_tn(:,:,:) * tmask(1:nlci,1:nlcj,1:jpk)
             tsn(1:nlci,1:nlcj,1:jpk,2) = temp_sn(:,:,:) * tmask(1:nlci,1:nlcj,1:jpk)
-            sshn(1:nlci,1:nlcj)        = temp_sshn(:,:) * tmask(1:nlci,1:nlcj,1)
+            ssh(1:nlci,1:nlcj,Nnn)     = temp_sshn(:,:) * tmask(1:nlci,1:nlcj,1)
 
-            ! Remove halo from tmask, tsn, sshn to prevent double obs counting
+            ! Remove halo from tmask, tsn, ssh to prevent double obs counting
             IF (jpi > nlci) THEN
                 tmask(nlci+1:,:,:) = 0
                 tsn(nlci+1:,:,:,1) = 0
                 tsn(nlci+1:,:,:,2) = 0
-                sshn(nlci+1:,:) = 0
+                ssh(nlci+1:,:,Nnn) = 0
             END IF
             IF (jpj > nlcj) THEN
@@ -153 +153 @@
                 tsn(:,nlcj+1:,:,1) = 0
                 tsn(:,nlcj+1:,:,2) = 0
-                sshn(:,nlcj+1:) = 0
+                ssh(:,nlcj+1:,Nnn) = 0
             END IF
 
@@ -161 +161 @@
             ! Mark all as missing data
             tsn(:,:,:,:) = fbrmdi
-            sshn(:,:) = fbrmdi
+            ssh(:,:,Nnn) = fbrmdi
          ENDIF
          ! Close netcdf file

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/TOP/oce_trc.F90

@@ -8 +8 @@
    !!----------------------------------------------------------------------
    !                                            !* Domain size *
-   USE par_oce , ONLY :   jpi      =>   jpi        !: first  dimension of grid --> i 
-   USE par_oce , ONLY :   jpj      =>   jpj        !: second dimension of grid --> j  
-   USE par_oce , ONLY :   jpk      =>   jpk        !: number of levels  
-   USE par_oce , ONLY :   jpim1    =>   jpim1      !: jpi - 1
-   USE par_oce , ONLY :   jpjm1    =>   jpjm1      !: jpj - 1 
-   USE par_oce , ONLY :   jpkm1    =>   jpkm1      !: jpk - 1  
-   USE par_oce , ONLY :   jpij     =>   jpij       !: jpi x jpj
-   USE par_oce , ONLY :   jp_tem   =>   jp_tem     !: indice for temperature
-   USE par_oce , ONLY :   jp_sal   =>   jp_sal     !: indice for salinity
+   USE par_oce, ONLY :   jpi      =>   jpi        !: first  dimension of grid --> i 
+   USE par_oce, ONLY :   jpj      =>   jpj        !: second dimension of grid --> j  
+   USE par_oce, ONLY :   jpk      =>   jpk        !: number of levels  
+   USE par_oce, ONLY :   jpim1    =>   jpim1      !: jpi - 1
+   USE par_oce, ONLY :   jpjm1    =>   jpjm1      !: jpj - 1 
+   USE par_oce, ONLY :   jpkm1    =>   jpkm1      !: jpk - 1  
+   USE par_oce, ONLY :   jpij     =>   jpij       !: jpi x jpj
+   USE par_oce, ONLY :   jp_tem   =>   jp_tem     !: indice for temperature
+   USE par_oce, ONLY :   jp_sal   =>   jp_sal     !: indice for salinity
+
+   USE par_oce, ONLY :   Nt       =>   Nt         !: number of time-levels
+   USE par_oce, ONLY :   Nbb      =>   Nbb        !: index for before time-level
+   USE par_oce, ONLY :   Nnn      =>   Nnn        !: index for now      -    -
+   USE par_oce, ONLY :   Naa      =>   Naa        !: index for after    -    -  
 
    USE in_out_manager                           !* IO manager *
@@ -42 +47 @@
    USE oce , ONLY :   rhd     =>    rhd     !: in situ density anomalie rhd=(rho-rho0)/rho0 (no units)
    USE oce , ONLY :   hdivn   =>    hdivn   !: horizontal divergence (1/s)
-   USE oce , ONLY :   sshn    =>    sshn    !: sea surface height at t-point [m]   
-   USE oce , ONLY :   sshb    =>    sshb    !: sea surface height at t-point [m]   
-   USE oce , ONLY :   ssha    =>    ssha    !: sea surface height at t-point [m]   
+   USE oce , ONLY :   ssh     =>    ssh     !: sea surface height at t-point [m]   
    USE oce , ONLY :   rab_n   =>    rab_n   !: local thermal/haline expension ratio at T-points
 
@@ -94 +97 @@
    USE zdfmxl , ONLY :   hmlpt       =>   hmlpt       !: mixed layer depth at t-points (m)
 
+   !!======================================================================
 END MODULE oce_trc

NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/TOP/trcsub.F90

@@ -104 +104 @@
            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
            !
-           sshb_hold  (:,:) = sshn  (:,:)
+           sshb_hold  (:,:) = ssh   (:,:,Nnn)
            emp_b_hold (:,:) = emp_b (:,:)
            !
@@ -136 +136 @@
          ENDIF 
          !
-         sshn_tm  (:,:)         = sshn_tm  (:,:)         + sshn  (:,:) 
+         sshn_tm  (:,:)         = sshn_tm  (:,:)         + ssh   (:,:,Nnn) 
          rnf_tm   (:,:)         = rnf_tm   (:,:)         + rnf   (:,:) 
          h_rnf_tm (:,:)         = h_rnf_tm (:,:)         + h_rnf (:,:) 
@@ -168 +168 @@
             ENDIF
          ENDIF 
-         sshn_temp  (:,:)        = sshn  (:,:)
-         sshb_temp  (:,:)        = sshb  (:,:)
-         ssha_temp  (:,:)        = ssha  (:,:)
+         sshn_temp  (:,:)        = ssh   (:,:,Nnn)
+         sshb_temp  (:,:)        = ssh   (:,:,Nbb)
+         ssha_temp  (:,:)        = ssh   (:,:,Naa)
          rnf_temp   (:,:)        = rnf   (:,:)
          h_rnf_temp (:,:)        = h_rnf (:,:)
@@ -206 +206 @@
             ENDIF
          ENDIF 
-         sshn_tm  (:,:)          = sshn_tm    (:,:)       + sshn  (:,:) 
-         rnf_tm   (:,:)          = rnf_tm     (:,:)       + rnf   (:,:) 
-         h_rnf_tm (:,:)          = h_rnf_tm   (:,:)       + h_rnf (:,:) 
+         sshn_tm  (:,:)          = sshn_tm    (:,:)       + ssh   (:,:,Nnn) 
+         rnf_tm   (:,:)          = rnf_tm     (:,:)       + rnf   (:,:)
+         h_rnf_tm (:,:)          = h_rnf_tm   (:,:)       + h_rnf (:,:)
          hmld_tm  (:,:)          = hmld_tm    (:,:)       + hmld  (:,:)
          fr_i_tm  (:,:)          = fr_i_tm    (:,:)       + fr_i  (:,:)
@@ -216 +216 @@
          wndm_tm  (:,:)          = wndm_tm    (:,:)       + wndm  (:,:)
          !
-         sshn     (:,:)          = sshn_tm    (:,:) * r1_ndttrcp1 
-         sshb     (:,:)          = sshb_hold  (:,:)
+         ssh      (:,:,Nnn)      = sshn_tm    (:,:) * r1_ndttrcp1 
+         ssh      (:,:,Nbb)      = sshb_hold  (:,:)
          rnf      (:,:)          = rnf_tm     (:,:) * r1_ndttrcp1 
          h_rnf    (:,:)          = h_rnf_tm   (:,:) * r1_ndttrcp1 
@@ -324 +324 @@
          vslp_tm (:,:,:)     = vslp (:,:,:)
       ENDIF
-      sshn_tm  (:,:) = sshn  (:,:) 
+      sshn_tm  (:,:) = ssh   (:,:,Nnn) 
       rnf_tm   (:,:) = rnf   (:,:) 
       h_rnf_tm (:,:) = h_rnf (:,:) 
@@ -377 +377 @@
          vslp  (:,:,:)=  vslp_temp  (:,:,:)
       ENDIF
-      sshn  (:,:)     =  sshn_temp  (:,:)
-      sshb  (:,:)     =  sshb_temp  (:,:)
-      ssha  (:,:)     =  ssha_temp  (:,:)
+      ssh   (:,:,Nnn) =  sshn_temp  (:,:)
+      ssh   (:,:,Nbb) =  sshb_temp  (:,:)
+      ssh   (:,:,Naa) =  ssha_temp  (:,:)
       rnf   (:,:)     =  rnf_temp   (:,:)
       h_rnf (:,:)     =  h_rnf_temp (:,:)
@@ -417 +417 @@
       ENDIF
       !
-      sshb_hold  (:,:) = sshn  (:,:)
+      sshb_hold  (:,:) = ssh   (:,:,Nnn)
       emp_b_hold (:,:) = emp   (:,:)
-      sshn_tm    (:,:) = sshn  (:,:) 
+      sshn_tm    (:,:) = ssh   (:,:,Nnn) 
       rnf_tm     (:,:) = rnf   (:,:) 
       h_rnf_tm   (:,:) = h_rnf (:,:) 
@@ -449 +449 @@
       !!                ***  ROUTINE trc_sub_ssh  ***
       !!                   
-      !! ** Purpose :   compute the after ssh (ssha), the now vertical velocity
+      !! ** Purpose :   compute the after ssh, the now vertical velocity
       !!              and update the now vertical coordinate (ln_linssh=F).
       !!
@@ -497 +497 @@
       ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
       z1_2rho0 = 0.5 * r1_rho0
-      ssha(:,:) = (  sshb(:,:) - rDt * ( z1_2rho0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
+      ssh(:,:,Naa) = (  ssh(:,:,Nbb) - rDt * ( z1_2rho0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
 
       IF( .NOT.ln_dynspg_ts ) THEN
@@ -506 +506 @@
 #endif
          IF( ln_bdy ) THEN
-            ssha(:,:) = ssha(:,:) * bdytmask(:,:)
-            CALL lbc_lnk( ssha, 'T', 1. ) 
+            ssh(:,:,Naa) = ssh(:,:,Naa) * bdytmask(:,:)
+            CALL lbc_lnk( ssh(:,:,Naa), 'T', 1. ) 
          ENDIF
       ENDIF