Changeset 12590 for NEMO/branches/2020

Timestamp:

2020-03-23T22:16:19+01:00 (4 years ago)

Author:

techene

Message:

all: add e3 substitute, OCE/DOM/domzgr_substitute.h90: correct a bug for e3f

Location:

NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE

Files:

• DOM/domzgr_substitute.h90 (modified) (1 diff)
• DYN/divhor.F90 (modified) (6 diffs)
• DYN/dynvor.F90 (modified) (22 diffs)
• DYN/sshwzv.F90 (modified) (19 diffs)
• FLO/floblk.F90 (modified) (24 diffs)
• SBC/sbcrnf.F90 (modified) (7 diffs)
• SBC/sbcssm.F90 (modified) (6 diffs)
• TRA/traadv_cen.F90 (modified) (9 diffs)
• TRA/traadv_fct.F90 (modified) (36 diffs)
• TRA/traadv_mus.F90 (modified) (10 diffs)
• TRA/traadv_qck.F90 (modified) (19 diffs)
• TRA/traadv_ubs.F90 (modified) (17 diffs)
• TRA/trabbc.F90 (modified) (6 diffs)
• TRA/trabbl.F90 (modified) (5 diffs)
• TRA/traisf.F90 (modified) (6 diffs)
• TRA/traldf_iso.F90 (modified) (17 diffs)
• TRA/traldf_lap_blp.F90 (modified) (10 diffs)
• TRA/traldf_triad.F90 (modified) (11 diffs)
• TRA/tramle.F90 (modified) (1 diff)
• TRA/tranpc.F90 (modified) (14 diffs)
• TRA/traqsr.F90 (modified) (19 diffs)
• TRA/trasbc.F90 (modified) (13 diffs)
• TRA/trazdf.F90 (modified) (9 diffs)

NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domzgr_substitute.h90

@@ -15 +15 @@
 #   define  e3u(i,j,k,t)   (e3u_0(i,j,k)*(1.+r3u(i,j,t)*umask(i,j,k)))
 #   define  e3v(i,j,k,t)   (e3v_0(i,j,k)*(1.+r3v(i,j,t)*vmask(i,j,k)))
-#   define  e3f(i,j,k,t)   (e3f_0(i,j,k)*(1.+r3f(i,j,t)*fmask(i,j,k)))
+#   define  e3f(i,j,k)     (e3f_0(i,j,k)*(1.+r3f(i,j)*fmask(i,j,k)))
 #   define  e3w(i,j,k,t)   (e3w_0(i,j,k)*(1.+r3t(i,j,t)))
 #   define  e3uw(i,j,k,t)  (e3uw_0(i,j,k)*(1.+r3u(i,j,t)))

NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/divhor.F90

@@ -21 +21 @@
    USE dom_oce         ! ocean space and time domain
    USE sbc_oce, ONLY : ln_rnf      ! river runoff
-   USE sbcrnf , ONLY : sbc_rnf_div ! river runoff 
+   USE sbcrnf , ONLY : sbc_rnf_div ! river runoff
    USE isf_oce, ONLY : ln_isf      ! ice shelf
    USE isfhdiv, ONLY : isf_hdiv    ! ice shelf
-#if defined key_asminc   
+#if defined key_asminc
    USE asminc          ! Assimilation increment
 #endif
@@ -40 +40 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$ 
+   !! $Id$
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
@@ -50 +51 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE div_hor  ***
-      !!                    
+      !!
       !! ** Purpose :   compute the horizontal divergence at now time-step
       !!
       !! ** Method  :   the now divergence is computed as :
       !!         hdiv = 1/(e1e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
-      !!      and correct with runoff inflow (div_rnf) and cross land flow (div_cla) 
+      !!      and correct with runoff inflow (div_rnf) and cross land flow (div_cla)
       !!
       !! ** Action  : - update hdiv, the now horizontal divergence
@@ -78 +79 @@
       DO_3D_00_00( 1, jpkm1 )
          hdiv(ji,jj,jk) = (  e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * uu(ji  ,jj,jk,Kmm)      &
-            &               - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm)      &
-            &               + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * vv(ji,jj  ,jk,Kmm)      &
-            &               - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm)  )   &
+            &              - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm)      &
+            &              + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * vv(ji,jj  ,jk,Kmm)      &
+            &              - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm)  )   &
             &            * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
       END_3D
@@ -95 +96 @@
       IF( ln_rnf )   CALL sbc_rnf_div( hdiv, Kmm )                     !==  runoffs    ==!   (update hdiv field)
       !
-#if defined key_asminc 
+#if defined key_asminc
       IF( ln_sshinc .AND. ln_asmiau )   CALL ssh_asm_div( kt, Kbb, Kmm, hdiv )   !==  SSH assimilation  ==!   (update hdiv field)
-      ! 
+      !
 #endif
       !
@@ -107 +108 @@
       !
    END SUBROUTINE div_hor
-   
+
    !!======================================================================
 END MODULE divhor

NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynvor.F90

@@ -15 +15 @@
    !!            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 
+   !!            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)
@@ -70 +70 @@
    INTEGER, PUBLIC, PARAMETER ::   np_MIX = 5   ! MIX scheme
 
-   INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity 
+   INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity
    !                               ! associated indices:
    INTEGER, PUBLIC, PARAMETER ::   np_COR = 1         ! Coriolis (planetary)
@@ -79 +79 @@
 
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2u_2        ! = di(e2u)/2          used in T-point metric term calculation
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1v_2        ! = dj(e1v)/2           -        -      -       - 
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1v_2        ! = dj(e1v)/2           -        -      -       -
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2v_2e1e2f   ! = di(e2u)/(2*e1e2f)  used in F-point metric term calculation
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1v)/(2*e1e2f)   -        -      -       - 
-   
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1v)/(2*e1e2f)   -        -      -       -
+
    REAL(wp) ::   r1_4  = 0.250_wp         ! =1/4
    REAL(wp) ::   r1_8  = 0.125_wp         ! =1/8
    REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! 1/12
-   
+
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
+
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -103 +105 @@
       !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend
       !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
-      !!               and planetary vorticity trends) and send them to trd_dyn 
+      !!               and planetary vorticity trends) and send them to trd_dyn
       !!               for futher diagnostics (l_trddyn=T)
       !!----------------------------------------------------------------------
@@ -193 +195 @@
       !!                  ***  ROUTINE vor_enT  ***
       !!
-      !! ** Purpose :   Compute the now total vorticity trend and add it to 
+      !! ** 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) 
+      !! ** Method  :   Trend evaluated using now fields (centered in time)
       !!       and t-point evaluation of vorticity (planetary and relative).
       !!       conserves the horizontal kinetic energy.
-      !!         The general trend of momentum is increased due to the vorticity 
+      !!         The general trend of momentum is increased due to the vorticity
       !!       term which is given by:
       !!          voru = 1/bu  mj[ ( mi(mj(bf*rvor))+bt*f_t)/e3t  mj[vn] ]
@@ -233 +235 @@
                   &             - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
             END_2D
-            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity 
+            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity
                DO_2D_10_10
                   zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
@@ -248 +250 @@
                   &              - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)   ) * r1_e1e2f(ji,jj)
             END_2D
-            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity 
+            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity
                DO_2D_10_10
                   zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
@@ -269 +271 @@
             DO_2D_01_01
                zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)   &
-                  &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
+                  &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) &
+                  &                  * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
             END_2D
          CASE ( np_MET )                           !* metric term
             DO_2D_01_01
-               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(ji,jj,jk,Kmm)
+               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(ji,jj,jk,Kmm)
             END_2D
          CASE ( np_CRV )                           !* Coriolis + relative vorticity
             DO_2D_01_01
-               zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)    &
-                  &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
+               zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)      &
+                  &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) &
+                  &                                 * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
             END_2D
          CASE ( np_CME )                           !* Coriolis + metric
             DO_2D_01_01
-               zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(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(ji,jj,jk,Kmm)
+               zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(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(ji,jj,jk,Kmm)
             END_2D
          CASE DEFAULT                                             ! error
@@ -298 +304 @@
                !
             pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)                    &
-               &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   & 
-               &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   ) 
+               &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   &
+               &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   )
          END_2D
          !                                             ! ===============
@@ -311 +317 @@
       !!                  ***  ROUTINE vor_ene  ***
       !!
-      !! ** Purpose :   Compute the now total vorticity trend and add it to 
+      !! ** 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) 
+      !! ** 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 
+      !!         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 pvv(:,:,:,Kmm)) ]
@@ -350 +356 @@
          SELECT CASE( kvor )                 !==  vorticity considered  ==!
          CASE ( np_COR )                           !* Coriolis (planetary vorticity)
-            zwz(:,:) = ff_f(:,:) 
+            zwz(:,:) = ff_f(:,:)
          CASE ( np_RVO )                           !* relative vorticity
             DO_2D_10_10
@@ -396 +402 @@
             zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
             pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
-            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
+            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
          END_2D
          !                                             ! ===============
@@ -446 +452 @@
          SELECT CASE( kvor )                 !==  vorticity considered  ==!
          CASE ( np_COR )                           !* Coriolis (planetary vorticity)
-            zwz(:,:) = ff_f(:,:) 
+            zwz(:,:) = ff_f(:,:)
          CASE ( np_RVO )                           !* relative vorticity
             DO_2D_10_10
@@ -504 +510 @@
       !!                ***  ROUTINE vor_een  ***
       !!
-      !! ** Purpose :   Compute the now total vorticity trend and add it to 
+      !! ** 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 
+      !! ** 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)
@@ -545 +551 @@
          CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
             DO_2D_10_10
-               ze3f = (  e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
-                  &    + e3t(ji,jj  ,jk,Kmm)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
+               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
+                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
+                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
+                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
                IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4._wp / ze3f
                ELSE                       ;   z1_e3f(ji,jj) = 0._wp
@@ -553 +561 @@
          CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
             DO_2D_10_10
-               ze3f = (  e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
-                  &    + e3t(ji,jj  ,jk,Kmm)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
+               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
+                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
+                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
+                  &    + e3t(ji+1,jj  ,jk,Kmm)*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)  )
@@ -644 +654 @@
       !!                ***  ROUTINE vor_eeT  ***
       !!
-      !! ** Purpose :   Compute the now total vorticity trend and add it to 
+      !! ** 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 
+      !! ** 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 
+      !!      The change consists in
       !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs).
       !!
@@ -667 +677 @@
       REAL(wp) ::   zua, zva       ! local scalars
       REAL(wp) ::   zmsk, z1_e3t   ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zwx , zwy 
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zwx , zwy
       REAL(wp), DIMENSION(jpi,jpj)     ::   ztnw, ztne, ztsw, ztse
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz
@@ -827 +837 @@
       !
       IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
-      !                      
+      !
       IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot)
       ncor = np_COR                       ! planetary vorticity
@@ -836 +846 @@
          ntot = np_COR        !     -         -
       CASE( np_VEC_c2  )
-         IF(lwp) WRITE(numout,*) '   ==>>>   vector form dynamics : total vorticity = Coriolis + relative vorticity' 
+         IF(lwp) WRITE(numout,*) '   ==>>>   vector form dynamics : total vorticity = Coriolis + relative vorticity'
          nrvm = np_RVO        ! relative vorticity
-         ntot = np_CRV        ! relative + planetary vorticity         
+         ntot = np_CRV        ! relative + planetary vorticity
       CASE( np_FLX_c2 , np_FLX_ubs  )
          IF(lwp) WRITE(numout,*) '   ==>>>   flux form dynamics : total vorticity = Coriolis + metric term'
@@ -863 +873 @@
          !
       END SELECT
-      
+
       IF(lwp) THEN                   ! Print the choice
          WRITE(numout,*)
@@ -873 +883 @@
          CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)'
          CASE( np_MIX )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)'
-         END SELECT         
+         END SELECT
       ENDIF
       !

NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/sshwzv.F90

@@ -1 @@
-MODULE sshwzv   
+MODULE sshwzv
    !!==============================================================================
    !!                       ***  MODULE  sshwzv  ***
@@ -5 +5 @@
    !!==============================================================================
    !! History :  3.1  !  2009-02  (G. Madec, M. Leclair)  Original code
-   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA 
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA
    !!             -   !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
    !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
@@ -20 +20 @@
    USE oce            ! ocean dynamics and tracers variables
    USE isf_oce        ! ice shelf
-   USE dom_oce        ! ocean space and time domain variables 
+   USE dom_oce        ! ocean space and time domain variables
    USE sbc_oce        ! surface boundary condition: ocean
    USE domvvl         ! Variable volume
@@ -31 +31 @@
 #endif
    !
-   USE iom 
+   USE iom
    USE in_out_manager ! I/O manager
    USE restart        ! only for lrst_oce
@@ -50 +50 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
+
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -60 +62 @@
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE ssh_nxt  ***
-      !!                   
+      !!
       !! ** Purpose :   compute the after ssh (ssh(Kaa))
       !!
@@ -74 +76 @@
       INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! time level index
       REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) ::   pssh           ! sea-surface height
-      ! 
+      !
       INTEGER  ::   jk            ! dummy loop indice
       REAL(wp) ::   z2dt, zcoef   ! local scalars
@@ -108 +110 @@
       ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
       ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
-      ! 
+      !
       pssh(:,:,Kaa) = (  pssh(:,:,Kbb) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
       !
@@ -131 +133 @@
    END SUBROUTINE ssh_nxt
 
-   
+
    SUBROUTINE wzv( kt, Kbb, Kmm, pww, Kaa )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE wzv  ***
-      !!                   
+      !!
       !! ** Purpose :   compute the now vertical velocity
       !!
-      !! ** Method  : - Using the incompressibility hypothesis, the vertical 
-      !!      velocity is computed by integrating the horizontal divergence  
+      !! ** Method  : - Using the incompressibility hypothesis, the vertical
+      !!      velocity is computed by integrating the horizontal divergence
       !!      from the bottom to the surface minus the scale factor evolution.
       !!        The boundary conditions are w=0 at the bottom (no flux) and.
@@ -172 +174 @@
       !
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
-         ALLOCATE( zhdiv(jpi,jpj,jpk) ) 
+         ALLOCATE( zhdiv(jpi,jpj,jpk) )
          !
          DO jk = 1, jpkm1
@@ -186 +188 @@
          DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
             ! computation of w
-            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk) + zhdiv(:,:,jk)    &
-               &                         + z1_2dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) )     ) * tmask(:,:,jk)
+            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)    &
+               &                         + zhdiv(:,:,jk)                       &
+               &                         + z1_2dt * ( e3t(:,:,jk,Kaa)          &
+               &                                    - e3t(:,:,jk,Kbb) )  ) * tmask(:,:,jk)
          END DO
          !          IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
-         DEALLOCATE( zhdiv ) 
+         DEALLOCATE( zhdiv )
       ELSE   ! z_star and linear free surface cases
          DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
             ! computation of w
-            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)                 &
-               &                         + z1_2dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) )  ) * tmask(:,:,jk)
+            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)    &
+               &                         + z1_2dt * ( e3t(:,:,jk,Kaa)          &
+               &                                    - e3t(:,:,jk,Kbb) )  ) * tmask(:,:,jk)
          END DO
       ENDIF
@@ -205 +210 @@
       ENDIF
       !
-#if defined key_agrif 
-      IF( .NOT. AGRIF_Root() ) THEN 
-         IF ((nbondi ==  1).OR.(nbondi == 2)) pww(nlci-1 , :     ,:) = 0.e0      ! east 
-         IF ((nbondi == -1).OR.(nbondi == 2)) pww(2      , :     ,:) = 0.e0      ! west 
-         IF ((nbondj ==  1).OR.(nbondj == 2)) pww(:      ,nlcj-1 ,:) = 0.e0      ! north 
-         IF ((nbondj == -1).OR.(nbondj == 2)) pww(:      ,2      ,:) = 0.e0      ! south 
-      ENDIF 
-#endif 
+#if defined key_agrif
+      IF( .NOT. AGRIF_Root() ) THEN
+         IF ((nbondi ==  1).OR.(nbondi == 2)) pww(nlci-1 , :     ,:) = 0.e0      ! east
+         IF ((nbondi == -1).OR.(nbondi == 2)) pww(2      , :     ,:) = 0.e0      ! west
+         IF ((nbondj ==  1).OR.(nbondj == 2)) pww(:      ,nlcj-1 ,:) = 0.e0      ! north
+         IF ((nbondj == -1).OR.(nbondj == 2)) pww(:      ,2      ,:) = 0.e0      ! south
+      ENDIF
+#endif
       !
       IF( ln_timing )   CALL timing_stop('wzv')
@@ -274 +279 @@
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE wAimp  ***
-      !!                   
+      !!
       !! ** Purpose :   compute the Courant number and partition vertical velocity
       !!                if a proportion needs to be treated implicitly
       !!
-      !! ** Method  : - 
+      !! ** Method  : -
       !!
       !! ** action  :   ww      : now vertical velocity (to be handled explicitly)
@@ -284 +289 @@
       !!
       !! Reference  : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent
-      !!              implicit scheme for vertical advection in oceanic modeling. 
+      !!              implicit scheme for vertical advection in oceanic modeling.
       !!              Ocean Modelling, 91, 38-69.
       !!----------------------------------------------------------------------
@@ -312 +317 @@
             z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
             ! 2*rdt and not r2dt (for restartability)
-            Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )                       &  
+            Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )                       &
                &                             + ( MAX( e2u(ji  ,jj)*e3u(ji  ,jj,jk,Kmm)*uu(ji  ,jj,jk,Kmm) + un_td(ji  ,jj,jk), 0._wp ) -   &
                &                                 MIN( e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm) + un_td(ji-1,jj,jk), 0._wp ) )   &
@@ -325 +330 @@
             z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
             ! 2*rdt and not r2dt (for restartability)
-            Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )   & 
+            Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )   &
                &                             + ( MAX( e2u(ji  ,jj)*e3u(ji  ,jj,jk,Kmm)*uu(ji  ,jj,jk,Kmm), 0._wp ) -   &
                &                                 MIN( e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm), 0._wp ) )   &
@@ -344 +349 @@
             zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) )
 ! alt:
-!                  IF ( ww(ji,jj,jk) > 0._wp ) THEN 
-!                     zCu =  Cu_adv(ji,jj,jk) 
+!                  IF ( ww(ji,jj,jk) > 0._wp ) THEN
+!                     zCu =  Cu_adv(ji,jj,jk)
 !                  ELSE
 !                     zCu =  Cu_adv(ji,jj,jk-1)
-!                  ENDIF 
+!                  ENDIF
             !
             IF( zCu <= Cu_min ) THEN              !<-- Fully explicit
@@ -365 +370 @@
             Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient below and in stp_ctl
          END_3D
-         Cu_adv(:,:,1) = 0._wp 
+         Cu_adv(:,:,1) = 0._wp
       ELSE
          ! Fully explicit everywhere
@@ -371 +376 @@
          wi    (:,:,:) = 0._wp
       ENDIF
-      CALL iom_put("wimp",wi) 
+      CALL iom_put("wimp",wi)
       CALL iom_put("wi_cff",Cu_adv)
       CALL iom_put("wexp",ww)

NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/FLO/floblk.F90

@@ -20 +20 @@
    PUBLIC   flo_blk    ! routine called by floats.F90
 
+#  include "domzgr_substitute.h90"
+
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$ 
+   !! $Id$
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
@@ -30 +32 @@
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE flo_blk  ***
-      !!           
+      !!
       !! ** Purpose :   Compute the geographical position,latitude, longitude
       !!      and depth of each float at each time step.
-      !! 
+      !!
       !! ** Method  :   The position of a float is computed with Bruno Blanke
       !!      algorithm. We need to know the velocity field, the old positions
@@ -47 +49 @@
          zuoutfl,zvoutfl,zwoutfl,   &     ! transport across the ouput face
          zvol,                      &     ! volume of the mesh
-         zsurfz,                    &     ! surface of the face of the mesh 
+         zsurfz,                    &     ! surface of the face of the mesh
          zind
 
@@ -53 +55 @@
 
       INTEGER  , DIMENSION ( jpnfl )  ::   iil, ijl, ikl                   ! index of nearest mesh
-      INTEGER  , DIMENSION ( jpnfl )  ::   iiloc , ijloc              
+      INTEGER  , DIMENSION ( jpnfl )  ::   iiloc , ijloc
       INTEGER  , DIMENSION ( jpnfl )  ::   iiinfl, ijinfl, ikinfl          ! index of input mesh of the float.
       INTEGER  , DIMENSION ( jpnfl )  ::   iioutfl, ijoutfl, ikoutfl       ! index of output mesh of the float.
-      REAL(wp) , DIMENSION ( jpnfl )  ::   zgifl, zgjfl, zgkfl             ! position of floats, index on 
+      REAL(wp) , DIMENSION ( jpnfl )  ::   zgifl, zgjfl, zgkfl             ! position of floats, index on
       !                                                                         ! velocity mesh.
       REAL(wp) , DIMENSION ( jpnfl )  ::    ztxfl, ztyfl, ztzfl            ! time for a float to quit the mesh
-      !                                                                         ! across one of the face x,y and z 
-      REAL(wp) , DIMENSION ( jpnfl )  ::    zttfl                          ! time for a float to quit the mesh 
-      REAL(wp) , DIMENSION ( jpnfl )  ::    zagefl                         ! time during which, trajectorie of 
+      !                                                                         ! across one of the face x,y and z
+      REAL(wp) , DIMENSION ( jpnfl )  ::    zttfl                          ! time for a float to quit the mesh
+      REAL(wp) , DIMENSION ( jpnfl )  ::    zagefl                         ! time during which, trajectorie of
       !                                                                         ! the float has been computed
-      REAL(wp) , DIMENSION ( jpnfl )  ::   zagenewfl                       ! new age of float after calculation 
+      REAL(wp) , DIMENSION ( jpnfl )  ::   zagenewfl                       ! new age of float after calculation
       !                                                                         ! of new position
       REAL(wp) , DIMENSION ( jpnfl )  ::   zufl, zvfl, zwfl                ! interpolated vel. at float position
@@ -77 +79 @@
 
       ! Initialisation of parameters
-      
+
       DO jfl = 1, jpnfl
          ! ages of floats are put at zero
          zagefl(jfl) = 0.
-         ! index on the velocity grid 
-         ! We considere k coordinate negative, with this transformation 
-         ! the computation in the 3 direction is the same. 
+         ! index on the velocity grid
+         ! We considere k coordinate negative, with this transformation
+         ! the computation in the 3 direction is the same.
          zgifl(jfl) = tpifl(jfl) - 0.5
          zgjfl(jfl) = tpjfl(jfl) - 0.5
          zgkfl(jfl) = MIN(-1.,-(tpkfl(jfl)))
-         ! surface drift every 10 days 
+         ! surface drift every 10 days
          IF( ln_argo ) THEN
             IF( MOD(kt,150) >= 146 .OR. MOD(kt,150) == 0 )  zgkfl(jfl) = -1.
@@ -96 +98 @@
          ikl(jfl) =     INT(zgkfl(jfl))
       END DO
-       
+
       iloop = 0
 222   DO jfl = 1, jpnfl
@@ -104 +106 @@
             iiloc(jfl) = iil(jfl) - mig(1) + 1
             ijloc(jfl) = ijl(jfl) - mjg(1) + 1
-# else 
+# else
             iiloc(jfl) = iil(jfl)
             ijloc(jfl) = ijl(jfl)
 # endif
-            
-            ! compute the transport across the mesh where the float is.            
-!!bug (gm) change e3t into e3. but never checked 
-            zsurfx(1) = e2u(iiloc(jfl)-1,ijloc(jfl)  ) * e3u(iiloc(jfl)-1,ijloc(jfl)  ,-ikl(jfl),Kmm)
-            zsurfx(2) = e2u(iiloc(jfl)  ,ijloc(jfl)  ) * e3u(iiloc(jfl)  ,ijloc(jfl)  ,-ikl(jfl),Kmm)
-            zsurfy(1) = e1v(iiloc(jfl)  ,ijloc(jfl)-1) * e3v(iiloc(jfl)  ,ijloc(jfl)-1,-ikl(jfl),Kmm)
-            zsurfy(2) = e1v(iiloc(jfl)  ,ijloc(jfl)  ) * e3v(iiloc(jfl)  ,ijloc(jfl)  ,-ikl(jfl),Kmm)
+
+            ! compute the transport across the mesh where the float is.
+!!bug (gm) change e3t into e3. but never checked
+            zsurfx(1) = e2u(iiloc(jfl)-1,ijloc(jfl)  )    &
+            &         * e3u(iiloc(jfl)-1,ijloc(jfl)  ,-ikl(jfl),Kmm)
+            zsurfx(2) = e2u(iiloc(jfl)  ,ijloc(jfl)  )    &
+            &         * e3u(iiloc(jfl)  ,ijloc(jfl)  ,-ikl(jfl),Kmm)
+            zsurfy(1) = e1v(iiloc(jfl)  ,ijloc(jfl)-1)    &
+            &         * e3v(iiloc(jfl)  ,ijloc(jfl)-1,-ikl(jfl),Kmm)
+            zsurfy(2) = e1v(iiloc(jfl)  ,ijloc(jfl)  )    &
+            &         * e3v(iiloc(jfl)  ,ijloc(jfl)  ,-ikl(jfl),Kmm)
 
             ! for a isobar float zsurfz is put to zero. The vertical velocity will be zero too.
@@ -129 +135 @@
             zwoutfl=-(wb(iiloc(jfl),ijloc(jfl),- ikl(jfl)   )   &
                &   +  ww(iiloc(jfl),ijloc(jfl),- ikl(jfl)   ) )/2. *  zsurfz*nisobfl(jfl)
-            
-            ! interpolation of velocity field on the float initial position            
+
+            ! interpolation of velocity field on the float initial position
             zufl(jfl)=  zuinfl  + ( zgifl(jfl) - float(iil(jfl)-1) ) * ( zuoutfl - zuinfl)
             zvfl(jfl)=  zvinfl  + ( zgjfl(jfl) - float(ijl(jfl)-1) ) * ( zvoutfl - zvinfl)
             zwfl(jfl)=  zwinfl  + ( zgkfl(jfl) - float(ikl(jfl)-1) ) * ( zwoutfl - zwinfl)
-            
+
             ! faces of input and output
             ! u-direction
@@ -147 +153 @@
                iiinfl (jfl) = iil(jfl) - 1
             ENDIF
-            ! v-direction       
+            ! v-direction
             IF( zvfl(jfl) < 0. ) THEN
                ijoutfl(jfl) = ijl(jfl) - 1.
@@ -169 +175 @@
                ikinfl (jfl) = ikl(jfl) - 1.
             ENDIF
-            
+
             ! compute the time to go out the mesh across a face
             ! u-direction
@@ -203 +209 @@
                ENDIF
             ENDIF
-            ! w-direction        
-            IF( nisobfl(jfl) == 1. ) THEN 
+            ! w-direction
+            IF( nisobfl(jfl) == 1. ) THEN
                zwdfl (jfl) = zwoutfl - zwinfl
                zgkdfl(jfl) = float(ikoutfl(jfl) - ikinfl(jfl))
@@ -221 +227 @@
                ENDIF
             ENDIF
-            
+
             ! the time to go leave the mesh is the smallest time
-                   
-            IF( nisobfl(jfl) == 1. ) THEN 
+
+            IF( nisobfl(jfl) == 1. ) THEN
                zttfl(jfl) = MIN(ztxfl(jfl),ztyfl(jfl),ztzfl(jfl))
             ELSE
@@ -231 +237 @@
             ! new age of the FLOAT
             zagenewfl(jfl) = zagefl(jfl) + zttfl(jfl)*zvol
-            ! test to know if the "age" of the float is not bigger than the 
+            ! test to know if the "age" of the float is not bigger than the
             ! time step
             IF( zagenewfl(jfl) > rdt ) THEN
@@ -237 +243 @@
                zagenewfl(jfl) = rdt
             ENDIF
-            
+
             ! In the "minimal" direction we compute the index of new mesh
             ! on i-direction
@@ -250 +256 @@
                iiinfl(jfl) = ind
             ELSE
-               IF( ABS(zudfl(jfl)) >= 1.E-5 ) THEN 
+               IF( ABS(zudfl(jfl)) >= 1.E-5 ) THEN
                   zgifl(jfl) = zgifl(jfl) + zgidfl(jfl)*zufl(jfl)    &
                      &       * ( EXP( zudfl(jfl)/zgidfl(jfl)*zttfl(jfl) ) - 1. ) /  zudfl(jfl)
@@ -268 +274 @@
                ijinfl(jfl) = ind
             ELSE
-               IF( ABS(zvdfl(jfl)) >= 1.E-5 ) THEN 
+               IF( ABS(zvdfl(jfl)) >= 1.E-5 ) THEN
                   zgjfl(jfl) = zgjfl(jfl)+zgjdfl(jfl)*zvfl(jfl)   &
                      &       * ( EXP(zvdfl(jfl)/zgjdfl(jfl)*zttfl(jfl)) - 1. ) /  zvdfl(jfl)
@@ -287 +293 @@
                   ikinfl(jfl) = ind
                ELSE
-                  IF( ABS(zwdfl(jfl)) >= 1.E-5 ) THEN 
+                  IF( ABS(zwdfl(jfl)) >= 1.E-5 ) THEN
                      zgkfl(jfl) = zgkfl(jfl)+zgkdfl(jfl)*zwfl(jfl)    &
                         &       * ( EXP(zwdfl(jfl)/zgkdfl(jfl)*zttfl(jfl)) - 1. ) /  zwdfl(jfl)
@@ -295 +301 @@
                ENDIF
             ENDIF
-            
+
             ! coordinate of the new point on the temperature grid
-            
+
             iil(jfl) = MAX(iiinfl(jfl),iioutfl(jfl))
             ijl(jfl) = MAX(ijinfl(jfl),ijoutfl(jfl))
@@ -306 +312 @@
 !!Alexcadm     .    ikoutfl(jfl),ztxfl(jfl),ztyfl(jfl)
 !!Alexcadm     .     ,ztzfl(jfl),zgifl(jfl),
-!!Alexcadm     .  zgjfl(jfl) 
+!!Alexcadm     .  zgjfl(jfl)
 !!Alexcadm  IF (jfl == 910) write(*,*)'Flotteur 910',
 !!Alexcadm     .    iiinfl(jfl),iioutfl(jfl),ijinfl(jfl)
@@ -312 +318 @@
 !!Alexcadm     .    ikoutfl(jfl),ztxfl(jfl),ztyfl(jfl)
 !!Alexcadm     .     ,ztzfl(jfl),zgifl(jfl),
-!!Alexcadm     .  zgjfl(jfl) 
+!!Alexcadm     .  zgjfl(jfl)
             ! reinitialisation of the age of FLOAT
             zagefl(jfl) = zagenewfl(jfl)
@@ -327 +333 @@
 # endif
       END DO
-      
+
       ! synchronisation
       CALL mpp_sum( 'floblk', zgifl , jpnfl )   ! sums over the global domain
@@ -335 +341 @@
       CALL mpp_sum( 'floblk', iil   , jpnfl )
       CALL mpp_sum( 'floblk', ijl   , jpnfl )
-      
+
       ! Test to know if a  float hasn't integrated enought time
       IF( ln_argo ) THEN
@@ -361 +367 @@
 !!Alexcadm     .       tpkfl(jpnfl),zufl(jpnfl),zvfl(jpnfl),zwfl(jpnfl)
       IF( ifin == 0 ) THEN
-         iloop = iloop + 1 
+         iloop = iloop + 1
          GO TO 222
       ENDIF
@@ -369 +375 @@
 
    !!======================================================================
-END MODULE floblk 
+END MODULE floblk

NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcrnf.F90

@@ -34 +34 @@
    PUBLIC   sbc_rnf_alloc ! called in sbcmod module
    PUBLIC   sbc_rnf_init  ! called in sbcmod module
-   
+
    !                                                !!* namsbc_rnf namelist *
    CHARACTER(len=100)         ::   cn_dir            !: Root directory for location of rnf files
@@ -58 +58 @@
    LOGICAL , PUBLIC ::   l_rnfcpl = .false.   !: runoffs recieved from oasis
    INTEGER , PUBLIC ::   nkrnf = 0            !: nb of levels over which Kz is increased at river mouths
-   
+
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rnfmsk              !: river mouth mask (hori.)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)     ::   rnfmsk_z            !: river mouth mask (vert.)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   h_rnf               !: depth of runoff in m
    INTEGER,  PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   nk_rnf              !: depth of runoff in model levels
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rnf_tsc_b, rnf_tsc  !: before and now T & S runoff contents   [K.m/s & PSU.m/s]   
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rnf_tsc_b, rnf_tsc  !: before and now T & S runoff contents   [K.m/s & PSU.m/s]
 
    TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_rnf       ! structure: river runoff (file information, fields read)
    TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_i_rnf     ! structure: iceberg flux (file information, fields read)
-   TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_s_rnf     ! structure: river runoff salinity (file information, fields read)  
-   TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_t_rnf     ! structure: river runoff temperature (file information, fields read)  
- 
+   TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_s_rnf     ! structure: river runoff salinity (file information, fields read)
+   TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_t_rnf     ! structure: river runoff temperature (file information, fields read)
+
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -227 +228 @@
       ELSE                       !==   runoff put only at the surface   ==!
          h_rnf (:,:)   = e3t (:,:,1,Kmm)        ! update h_rnf to be depth of top box
-         phdivn(:,:,1) = phdivn(:,:,1) - ( rnf(:,:) + rnf_b(:,:) ) * zfact * r1_rau0 / e3t(:,:,1,Kmm)
+         phdivn(:,:,1) = phdivn(:,:,1) - ( rnf(:,:)+rnf_b(:,:) ) * zfact * r1_rau0 / e3t(:,:,1,Kmm)
       ENDIF
       !
@@ -249 +250 @@
       INTEGER           ::   ios           ! Local integer output status for namelist read
       INTEGER           ::   nbrec         ! temporary integer
-      REAL(wp)          ::   zacoef  
-      REAL(wp), DIMENSION(jpi,jpj,2) :: zrnfcl    
+      REAL(wp)          ::   zacoef
+      REAL(wp), DIMENSION(jpi,jpj,2) :: zrnfcl
       !!
       NAMELIST/namsbc_rnf/ cn_dir            , ln_rnf_depth, ln_rnf_tem, ln_rnf_sal, ln_rnf_icb,   &
@@ -261 +262 @@
       IF( sbc_rnf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_rnf_alloc : unable to allocate arrays' )
       !
-      IF( .NOT. ln_rnf ) THEN                      ! no specific treatment in vicinity of river mouths 
+      IF( .NOT. ln_rnf ) THEN                      ! no specific treatment in vicinity of river mouths
          ln_rnf_mouth  = .FALSE.                   ! default definition needed for example by sbc_ssr or by tra_adv_muscl
          nkrnf         = 0
@@ -297 +298 @@
       !                                   ! ==================
       !
-      IF( .NOT. l_rnfcpl ) THEN                    
+      IF( .NOT. l_rnfcpl ) THEN
          ALLOCATE( sf_rnf(1), STAT=ierror )         ! Create sf_rnf structure (runoff inflow)
          IF(lwp) WRITE(numout,*)
@@ -352 +353 @@
          IF(lwp) WRITE(numout,*) '   ==>>>   runoffs depth read in a file'
          rn_dep_file = TRIM( cn_dir )//TRIM( sn_dep_rnf%clname )
-         IF( .NOT. sn_dep_rnf%ln_clim ) THEN   ;   WRITE(rn_dep_file, '(a,"_y",i4)' ) TRIM( rn_dep_file ), nyear    ! add year 
-            IF( sn_dep_rnf%cltype == 'monthly' )   WRITE(rn_dep_file, '(a,"m",i2)'  ) TRIM( rn_dep_file ), nmonth   ! add month 
+         IF( .NOT. sn_dep_rnf%ln_clim ) THEN   ;   WRITE(rn_dep_file, '(a,"_y",i4)' ) TRIM( rn_dep_file ), nyear    ! add year
+            IF( sn_dep_rnf%cltype == 'monthly' )   WRITE(rn_dep_file, '(a,"m",i2)'  ) TRIM( rn_dep_file ), nmonth   ! add month
          ENDIF
          CALL iom_open ( rn_dep_file, inum )                           ! open file

NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcssm.F90

@@ -10 +10 @@
 
    !!----------------------------------------------------------------------
-   !!   sbc_ssm       : calculate sea surface mean currents, temperature,  
+   !!   sbc_ssm       : calculate sea surface mean currents, temperature,
    !!                   and salinity over nn_fsbc time-step
    !!----------------------------------------------------------------------
@@ -31 +31 @@
 
    LOGICAL, SAVE ::   l_ssm_mean = .FALSE.   ! keep track of whether means have been read from restart file
-   
+
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -42 +43 @@
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE sbc_oce  ***
-      !!                     
+      !!
       !! ** Purpose :   provide ocean surface variable to sea-surface boundary
-      !!                condition computation 
-      !!                
-      !! ** Method  :   compute mean surface velocity (2 components at U and 
+      !!                condition computation
+      !!
+      !! ** Method  :   compute mean surface velocity (2 components at U and
       !!      V-points) [m/s], temperature [Celsius] and salinity [psu] over
       !!      the periode (kt - nn_fsbc) to kt
@@ -200 +201 @@
          !
       ELSE
-         !               
+         !
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'sbc_ssm_init : sea surface mean fields'
@@ -222 +223 @@
             !
             IF( zf_sbc /= REAL( nn_fsbc, wp ) ) THEN      ! nn_fsbc has changed between 2 runs
-               IF(lwp) WRITE(numout,*) '   restart with a change in the frequency of mean from ', zf_sbc, ' to ', nn_fsbc 
-               zcoef = REAL( nn_fsbc - 1, wp ) / zf_sbc 
-               ssu_m(:,:) = zcoef * ssu_m(:,:) 
+               IF(lwp) WRITE(numout,*) '   restart with a change in the frequency of mean from ', zf_sbc, ' to ', nn_fsbc
+               zcoef = REAL( nn_fsbc - 1, wp ) / zf_sbc
+               ssu_m(:,:) = zcoef * ssu_m(:,:)
                ssv_m(:,:) = zcoef * ssv_m(:,:)
                sst_m(:,:) = zcoef * sst_m(:,:)
@@ -252 +253 @@
       ENDIF
       !
-      IF( .NOT. ln_traqsr )   fraqsr_1lev(:,:) = 1._wp   ! default definition: qsr 100% in the fisrt level 
+      IF( .NOT. ln_traqsr )   fraqsr_1lev(:,:) = 1._wp   ! default definition: qsr 100% in the fisrt level
       !
       IF( lwxios.AND.nn_fsbc > 1 ) THEN

NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traadv_cen.F90

@@ -13 +13 @@
    USE dom_oce        ! ocean space and time domain
    USE eosbn2         ! equation of state
-   USE traadv_fct     ! acces to routine interp_4th_cpt 
+   USE traadv_fct     ! acces to routine interp_4th_cpt
    USE trd_oce        ! trends: ocean variables
-   USE trdtra         ! trends manager: tracers 
+   USE trdtra         ! trends manager: tracers
    USE diaptr         ! poleward transport diagnostics
    USE diaar5         ! AR5 diagnostics
@@ -28 +28 @@
 
    PUBLIC   tra_adv_cen   ! called by traadv.F90
-   
+
    REAL(wp) ::   r1_6 = 1._wp / 6._wp   ! =1/6
 
@@ -37 +37 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -45 +46 @@
 
    SUBROUTINE tra_adv_cen( kt, kit000, cdtype, pU, pV, pW,     &
-      &                    Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v ) 
+      &                    Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv_cen  ***
-      !!                 
+      !!
       !! ** Purpose :   Compute the now trend due to the advection of tracers
       !!      and add it to the general trend of passive tracer equations.
       !!
       !! ** Method  :   The advection is evaluated by a 2nd or 4th order scheme
-      !!               using now fields (leap-frog scheme). 
+      !!               using now fields (leap-frog scheme).
       !!       kn_cen_h = 2  ==>> 2nd order centered scheme on the horizontal
       !!                = 4  ==>> 4th order    -        -       -      -
@@ -90 +91 @@
       l_ptr = .FALSE.
       IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )       l_trd = .TRUE.
-      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )    l_ptr = .TRUE. 
+      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )    l_ptr = .TRUE.
       IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
          &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE.
       !
-      !                    
+      !
       zwz(:,:, 1 ) = 0._wp       ! surface & bottom vertical flux set to zero for all tracers
       zwz(:,:,jpk) = 0._wp
@@ -150 +151 @@
             IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
                DO_2D_11_11
-                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm) 
+                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm)
                END_2D
             ELSE                                   ! no ice-shelf cavities (only ocean surface)
@@ -156 +157 @@
             ENDIF
          ENDIF
-         !               
+         !
          DO_3D_00_00( 1, jpkm1 )
             pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs)    &
                &             - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )    &
                &                + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )    &
-               &                + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)  ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+               &                + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)  ) &
+               &                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
          END_3D
          !                             ! trend diagnostics
@@ -169 +171 @@
             CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) )
          END IF
-         !                                 ! "Poleward" heat and salt transports 
+         !                                 ! "Poleward" heat and salt transports
          IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
          !                                 !  heat and salt transport
@@ -177 +179 @@
       !
    END SUBROUTINE tra_adv_cen
-   
+
    !!======================================================================
 END MODULE traadv_cen

NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traadv_fct.F90

@@ -10 +10 @@
    !!  tra_adv_fct    : update the tracer trend with a 3D advective trends using a 2nd or 4th order FCT scheme
    !!                   with sub-time-stepping in the vertical direction
-   !!  nonosc         : compute monotonic tracer fluxes by a non-oscillatory algorithm 
+   !!  nonosc         : compute monotonic tracer fluxes by a non-oscillatory algorithm
    !!  interp_4th_cpt : 4th order compact scheme for the vertical component of the advection
    !!----------------------------------------------------------------------
@@ -24 +24 @@
    !
    USE in_out_manager ! I/O manager
-   USE iom            ! 
+   USE iom            !
    USE lib_mpp        ! MPP library
-   USE lbclnk         ! ocean lateral boundary condition (or mpp link) 
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
 
    IMPLICIT NONE
@@ -46 +46 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -57 +58 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv_fct  ***
-      !! 
+      !!
       !! **  Purpose :   Compute the now trend due to total advection of tracers
       !!               and add it to the general trend of tracer equations
@@ -63 +64 @@
       !! **  Method  : - 2nd or 4th FCT scheme on the horizontal direction
       !!               (choice through the value of kn_fct)
-      !!               - on the vertical the 4th order is a compact scheme 
-      !!               - corrected flux (monotonic correction) 
+      !!               - on the vertical the 4th order is a compact scheme
+      !!               - corrected flux (monotonic correction)
       !!
       !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends
@@ -81 +82 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
       !
-      INTEGER  ::   ji, jj, jk, jn                           ! dummy loop indices  
+      INTEGER  ::   ji, jj, jk, jn                           ! dummy loop indices
       REAL(wp) ::   ztra                                     ! local scalar
       REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk, zC2t_u, zC4t_u   !   -      -
@@ -102 +103 @@
       ll_zAimp = .FALSE.
       IF( ( cdtype == 'TRA' .AND. l_trdtra  ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
-      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) )    l_ptr = .TRUE. 
+      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) )    l_ptr = .TRUE.
       IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR.  &
          &                         iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE.
@@ -111 +112 @@
       ENDIF
       !
-      IF( l_ptr ) THEN  
+      IF( l_ptr ) THEN
          ALLOCATE( zptry(jpi,jpj,jpk) )
          zptry(:,:,:) = 0._wp
       ENDIF
       !                          ! surface & bottom value : flux set to zero one for all
-      zwz(:,:, 1 ) = 0._wp            
+      zwz(:,:, 1 ) = 0._wp
       zwx(:,:,jpk) = 0._wp   ;   zwy(:,:,jpk) = 0._wp    ;    zwz(:,:,jpk) = 0._wp
       !
-      zwi(:,:,:) = 0._wp        
+      zwi(:,:,:) = 0._wp
       !
       ! If adaptive vertical advection, check if it is needed on this PE at this time
@@ -129 +130 @@
          ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk))
          DO_3D_00_00( 1, jpkm1 )
-            zwdia(ji,jj,jk) =  1._wp + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) / e3t(ji,jj,jk,Krhs)
+            zwdia(ji,jj,jk) =  1._wp + p2dt * ( MAX( wi(ji,jj,jk) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) )   &
+            &                               / e3t(ji,jj,jk,Krhs)
             zwinf(ji,jj,jk) =  p2dt * MIN( wi(ji,jj,jk  ) , 0._wp ) / e3t(ji,jj,jk,Krhs)
             zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t(ji,jj,jk,Krhs)
@@ -138 +140 @@
          !
          !        !==  upstream advection with initial mass fluxes & intermediate update  ==!
-         !                    !* upstream tracer flux in the i and j direction 
+         !                    !* upstream tracer flux in the i and j direction
          DO_3D_10_10( 1, jpkm1 )
             ! upstream scheme
@@ -157 +159 @@
             IF( ln_isfcav ) THEN             ! top of the ice-shelf cavities and at the ocean surface
                DO_2D_11_11
-                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface 
+                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface
                END_2D
             ELSE                             ! no cavities: only at the ocean surface
@@ -163 +165 @@
             ENDIF
          ENDIF
-         !               
+         !
          DO_3D_00_00( 1, jpkm1 )
             !                             ! total intermediate advective trends
@@ -170 +172 @@
                &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj)
             !                             ! update and guess with monotonic sheme
-            pt(ji,jj,jk,jn,Krhs) =                     pt(ji,jj,jk,jn,Krhs) +        ztra   / e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
-            zwi(ji,jj,jk)    = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
+            pt(ji,jj,jk,jn,Krhs) =                   pt(ji,jj,jk,jn,Krhs) +       ztra   &
+               &                                  / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk)
+            zwi(ji,jj,jk)    = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) &
+               &                                  / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
          END_3D
-         
+
          IF ( ll_zAimp ) THEN
             CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 )
@@ -186 +190 @@
             DO_3D_00_00( 1, jpkm1 )
                pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) &
-                  &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+                  &                                     * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
             END_3D
             !
          END IF
-         !                
+         !
          IF( l_trd .OR. l_hst )  THEN             ! trend diagnostics (contribution of upstream fluxes)
             ztrdx(:,:,:) = zwx(:,:,:)   ;   ztrdy(:,:,:) = zwy(:,:,:)   ;   ztrdz(:,:,:) = zwz(:,:,:)
          END IF
          !                             ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( l_ptr )   zptry(:,:,:) = zwy(:,:,:) 
+         IF( l_ptr )   zptry(:,:,:) = zwy(:,:,:)
          !
          !        !==  anti-diffusive flux : high order minus low order  ==!
@@ -225 +229 @@
                zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! 2 x C2 interpolation of T at u- & v-points
                zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm)
-               !                                                  ! C4 minus upstream advective fluxes 
+               !                                                  ! C4 minus upstream advective fluxes
                zwx(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk)
                zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk)
@@ -245 +249 @@
                zC4t_u =  zC2t_u + r1_6 * ( ztu(ji-1,jj  ,jk) - ztu(ji+1,jj  ,jk) )
                zC4t_v =  zC2t_v + r1_6 * ( ztv(ji  ,jj-1,jk) - ztv(ji  ,jj+1,jk) )
-               !                                                  ! C4 minus upstream advective fluxes 
+               !                                                  ! C4 minus upstream advective fluxes
                zwx(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk)
                zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk)
@@ -251 +255 @@
             !
          END SELECT
-         !                      
+         !
          SELECT CASE( kn_fct_v )    !* vertical anti-diffusive fluxes (w-masked interior values)
          !
@@ -270 +274 @@
             zwz(:,:,1) = 0._wp   ! only ocean surface as interior zwz values have been w-masked
          ENDIF
-         !         
+         !
          IF ( ll_zAimp ) THEN
             DO_3D_00_00( 1, jpkm1 )
@@ -277 +281 @@
                   &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
                   &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj)
-               ztw(ji,jj,jk)  = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
+               ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs)*tmask(ji,jj,jk)
             END_3D
             !
@@ -316 +320 @@
             DO_3D_00_00( 1, jpkm1 )
                pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) &
-                  &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
-            END_3D
-         END IF         
+                  &                                     * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+            END_3D
+         END IF
          !
          IF( l_trd .OR. l_hst ) THEN   ! trend diagnostics // heat/salt transport
-            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< add anti-diffusive fluxes 
+            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< add anti-diffusive fluxes
             ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  !     to upstream fluxes
             ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  !
@@ -344 +348 @@
          DEALLOCATE( zwdia, zwinf, zwsup )
       ENDIF
-      IF( l_trd .OR. l_hst ) THEN 
+      IF( l_trd .OR. l_hst ) THEN
          DEALLOCATE( ztrdx, ztrdy, ztrdz )
       ENDIF
-      IF( l_ptr ) THEN 
+      IF( l_ptr ) THEN
          DEALLOCATE( zptry )
       ENDIF
@@ -357 +361 @@
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE nonosc  ***
-      !!     
-      !! **  Purpose :   compute monotonic tracer fluxes from the upstream 
-      !!       scheme and the before field by a nonoscillatory algorithm 
+      !!
+      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
+      !!       scheme and the before field by a nonoscillatory algorithm
       !!
       !! **  Method  :   ... ???
@@ -367 +371 @@
       !!       in-space based differencing for fluid
       !!----------------------------------------------------------------------
-      INTEGER                          , INTENT(in   ) ::   Kmm             ! time level index 
+      INTEGER                          , INTENT(in   ) ::   Kmm             ! time level index
       REAL(wp)                         , INTENT(in   ) ::   p2dt            ! tracer time-step
       REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   pbef, paft      ! before & after field
@@ -453 +457 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE interp_4th_cpt_org  ***
-      !! 
+      !!
       !! **  Purpose :   Compute the interpolation of tracer at w-point
       !!
@@ -464 +468 @@
       REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt
       !!----------------------------------------------------------------------
-      
+
       DO_3D_11_11( 3, jpkm1 )
          zwd (ji,jj,jk) = 4._wp
@@ -475 +479 @@
             zwi (ji,jj,jk) = 0._wp
             zws (ji,jj,jk) = 0._wp
-            zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )    
+            zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
          ENDIF
       END_3D
@@ -499 +503 @@
       END_2D
       DO_3D_11_11( 3, jpkm1 )
-         pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)             
+         pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)
       END_3D
 
@@ -508 +512 @@
          pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
       END_3D
-      !    
+      !
    END SUBROUTINE interp_4th_cpt_org
-   
+
 
    SUBROUTINE interp_4th_cpt( pt_in, pt_out )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE interp_4th_cpt  ***
-      !! 
+      !!
       !! **  Purpose :   Compute the interpolation of tracer at w-point
       !!
@@ -543 +547 @@
 !      CASE( np_CEN2 )   ! 2nd order centered  at top & bottom
 !      END SELECT
-!!gm  
+!!gm
       !
       IF ( ln_isfcav ) THEN            ! set level two values which may not be set in ISF case
@@ -561 +565 @@
          zwi (ji,jj,ikb) = 0._wp
          zws (ji,jj,ikb) = 0._wp
-         zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,ikb-1) + pt_in(ji,jj,ikb) )            
+         zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,ikb-1) + pt_in(ji,jj,ikb) )
       END_2D
       !
@@ -577 +581 @@
       END_2D
       DO_3D_00_00( 3, jpkm1 )
-         pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)             
+         pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)
       END_3D
 
@@ -586 +590 @@
          pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
       END_3D
-      !    
+      !
    END SUBROUTINE interp_4th_cpt
 
@@ -593 +597 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tridia_solver  ***
-      !! 
+      !!
       !! **  Purpose :   solve a symmetric 3diagonal system
       !!
       !! **  Method  :   solve M.t_out = RHS(t)  where M is a tri diagonal matrix ( jpk*jpk )
-      !!     
+      !!
       !!             ( D_1 U_1  0   0   0  )( t_1 )   ( RHS_1 )
       !!             ( L_2 D_2 U_2  0   0  )( t_2 )   ( RHS_2 )
@@ -603 +607 @@
       !!             (        ...          )( ... )   ( ...  )
       !!             (  0   0   0  L_k D_k )( t_k )   ( RHS_k )
-      !!     
+      !!
       !!        M is decomposed in the product of an upper and lower triangular matrix.
-      !!        The tri-diagonals matrix is given as input 3D arrays:   pD, pU, pL 
+      !!        The tri-diagonals matrix is given as input 3D arrays:   pD, pU, pL
       !!        (i.e. the Diagonal, the Upper diagonal, and the Lower diagonal).
       !!        The solution is pta.
@@ -613 +617 @@
       REAL(wp),DIMENSION(:,:,:), INTENT(in   ) ::   pRHS          ! Right-Hand-Side
       REAL(wp),DIMENSION(:,:,:), INTENT(  out) ::   pt_out        !!gm field at level=F(klev)
-      INTEGER                  , INTENT(in   ) ::   klev          ! =1 pt_out at w-level 
+      INTEGER                  , INTENT(in   ) ::   klev          ! =1 pt_out at w-level
       !                                                           ! =0 pt at t-level
       INTEGER ::   ji, jj, jk   ! dummy loop integers
@@ -633 +637 @@
       END_2D
       DO_3D_00_00( kstart+1, jpkm1 )
-         pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)             
+         pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)
       END_3D
 

NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traadv_mus.F90

@@ -29 +29 @@
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! distribued memory computing
-   USE lbclnk         ! ocean lateral boundary condition (or mpp link) 
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
 
    IMPLICIT NONE
@@ -36 +36 @@
 
    PUBLIC   tra_adv_mus   ! routine called by traadv.F90
-   
+
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   upsmsk   !: mixed upstream/centered scheme near some straits
    !                                                           !  and in closed seas (orca 2 and 1 configurations)
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xind     !: mixed upstream/centered index
-   
+
    LOGICAL  ::   l_trd   ! flag to compute trends
    LOGICAL  ::   l_ptr   ! flag to compute poleward transport
@@ -47 +47 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$ 
+   !! $Id$
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
@@ -64 +65 @@
       !!
       !! ** Method  : MUSCL scheme plus centered scheme at ocean boundaries
-      !!              ld_msc_ups=T : 
+      !!              ld_msc_ups=T :
       !!
       !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends
@@ -88 +89 @@
       REAL(wp) ::   zv, z0v, zzwy, z0w           !   -      -
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zslpx   ! 3D workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwy, zslpy   ! -      - 
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwy, zslpy   ! -      -
       !!----------------------------------------------------------------------
       !
@@ -112 +113 @@
                   &                 upsmsk(:,:)                ) * tmask(:,:,jk)   ! =>0 in some user defined area
             END DO
-         ENDIF 
-         !
-      ENDIF 
-      !      
+         ENDIF
+         !
+      ENDIF
+      !
       l_trd = .FALSE.
       l_hst = .FALSE.
       l_ptr = .FALSE.
       IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
-      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )   l_ptr = .TRUE. 
+      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )   l_ptr = .TRUE.
       IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
          &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE.
@@ -130 +131 @@
          !                                !-- first guess of the slopes
          zwx(:,:,jpk) = 0._wp                   ! bottom values
-         zwy(:,:,jpk) = 0._wp  
+         zwy(:,:,jpk) = 0._wp
          DO_3D_10_10( 1, jpkm1 )
             zwx(ji,jj,jk) = umask(ji,jj,jk) * ( pt(ji+1,jj,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
@@ -176 +177 @@
          DO_3D_00_00( 1, jpkm1 )
             pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )       &
-            &                                     + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )     &
-            &                                   * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+            &                                           +   zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )     &
+            &                           * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
          END_3D
          !                                ! trend diagnostics
@@ -184 +185 @@
             CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kbb) )
          END IF
-         !                                 ! "Poleward" heat and salt transports 
+         !                                 ! "Poleward" heat and salt transports
          IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
          !                                 !  heat transport
@@ -227 +228 @@
          !
          DO_3D_00_00( 1, jpkm1 )
-            pt(ji,jj,jk,jn,Krhs) =  pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+            pt(ji,jj,jk,jn,Krhs) =  pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )   &
+               &                                      * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
          END_3D
          !                                ! send trends for diagnostic

NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traadv_qck.F90

@@ -19 +19 @@
    USE trc_oce         ! share passive tracers/Ocean variables
    USE trd_oce         ! trends: ocean variables
-   USE trdtra          ! trends manager: tracers 
+   USE trdtra          ! trends manager: tracers
    USE diaptr          ! poleward transport diagnostics
    USE iom
@@ -26 +26 @@
    USE lib_mpp         ! distribued memory computing
    USE lbclnk          ! ocean lateral boundary condition (or mpp link)
-   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
 
    IMPLICIT NONE
@@ -41 +41 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -104 +105 @@
       l_ptr = .FALSE.
       IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE.
-      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE. 
+      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.
       !
       !
       !        ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
-      CALL tra_adv_qck_i( kt, cdtype, p2dt, pU, Kbb, Kmm, pt, kjpt, Krhs ) 
-      CALL tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs ) 
+      CALL tra_adv_qck_i( kt, cdtype, p2dt, pU, Kbb, Kmm, pt, kjpt, Krhs )
+      CALL tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs )
 
       !        ! vertical fluxes are computed with the 2nd order centered scheme
@@ -137 +138 @@
       DO jn = 1, kjpt                                            ! tracer loop
          !                                                       ! ===========
-         zfu(:,:,:) = 0._wp     ;   zfc(:,:,:) = 0._wp 
-         zfd(:,:,:) = 0._wp     ;   zwx(:,:,:) = 0._wp   
+         zfu(:,:,:) = 0._wp     ;   zfc(:,:,:) = 0._wp
+         zfd(:,:,:) = 0._wp     ;   zwx(:,:,:) = 0._wp
          !
 !!gm why not using a SHIFT instruction...
@@ -145 +146 @@
             zfd(ji,jj,jk) = pt(ji+1,jj,jk,jn,Kbb)        ! Downstream in the x-direction for the tracer
          END_3D
-         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions 
-         
+         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
+
          !
          ! Horizontal advective fluxes
          ! ---------------------------
          DO_3D_00_00( 1, jpkm1 )
-            zdir = 0.5 + SIGN( 0.5, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
-            zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T 
-         END_3D
-         !
-         DO_3D_00_00( 1, jpkm1 )
-            zdir = 0.5 + SIGN( 0.5, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
+            zdir = 0.5 + SIGN( 0.5, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
+            zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T
+         END_3D
+         !
+         DO_3D_00_00( 1, jpkm1 )
+            zdir = 0.5 + SIGN( 0.5, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
             zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
             zwx(ji,jj,jk)  = ABS( pU(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
@@ -162 +163 @@
             zfd(ji,jj,jk)  = zdir * pt(ji+1,jj,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji  ,jj,jk,jn,Kbb)  ! FD in the x-direction for T
          END_3D
-         !--- Lateral boundary conditions 
+         !--- Lateral boundary conditions
          CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1., zfc(:,:,:), 'T', 1.,  zwx(:,:,:), 'T', 1. )
 
@@ -172 +173 @@
             zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
          END_3D
-         CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1. )      ! Lateral boundary conditions 
+         CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1. )      ! Lateral boundary conditions
 
          !
          ! Tracer flux on the x-direction
-         DO jk = 1, jpkm1  
+         DO jk = 1, jpkm1
             !
             DO_2D_00_00
-               zdir = 0.5 + SIGN( 0.5, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
+               zdir = 0.5 + SIGN( 0.5, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
                !--- If the second ustream point is a land point
                !--- the flux is computed by the 1st order UPWIND scheme
@@ -226 +227 @@
       DO jn = 1, kjpt                                            ! tracer loop
          !                                                       ! ===========
-         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0  
-         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0     
-         !                                                  
-         DO jk = 1, jpkm1                                
-            !                                             
+         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0
+         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0
+         !
+         DO jk = 1, jpkm1
+            !
             !--- Computation of the ustream and downstream value of the tracer and the mask
             DO_2D_00_00
@@ -239 +240 @@
             END_2D
          END DO
-         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions 
-
-         
+         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
+
+
          !
          ! Horizontal advective fluxes
@@ -247 +248 @@
          !
          DO_3D_00_00( 1, jpkm1 )
-            zdir = 0.5 + SIGN( 0.5, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
-            zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T 
-         END_3D
-         !
-         DO_3D_00_00( 1, jpkm1 )
-            zdir = 0.5 + SIGN( 0.5, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
-            zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
+            zdir = 0.5 + SIGN( 0.5, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
+            zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
+         END_3D
+         !
+         DO_3D_00_00( 1, jpkm1 )
+            zdir = 0.5 + SIGN( 0.5, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
+            zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) )   &
+               &         * e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
             zwy(ji,jj,jk)  = ABS( pV(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
             zfc(ji,jj,jk)  = zdir * pt(ji,jj  ,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji,jj+1,jk,jn,Kbb)  ! FC in the x-direction for T
@@ -259 +261 @@
          END_3D
 
-         !--- Lateral boundary conditions 
+         !--- Lateral boundary conditions
          CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1., zfc(:,:,:), 'T', 1., zwy(:,:,:), 'T', 1. )
 
@@ -269 +271 @@
             zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
          END_3D
-         CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1. )    !--- Lateral boundary conditions 
+         CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1. )    !--- Lateral boundary conditions
          !
          ! Tracer flux on the x-direction
-         DO jk = 1, jpkm1  
+         DO jk = 1, jpkm1
             !
             DO_2D_00_00
-               zdir = 0.5 + SIGN( 0.5, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
+               zdir = 0.5 + SIGN( 0.5, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
                !--- If the second ustream point is a land point
                !--- the flux is computed by the 1st order UPWIND scheme
@@ -312 +314 @@
       CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
       INTEGER                                  , INTENT(in   ) ::   kjpt     ! number of tracers
-      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pW      ! vertical velocity 
+      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pW      ! vertical velocity
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
       !
@@ -332 +334 @@
             IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
                DO_2D_11_11
-                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm)   ! linear free surface 
+                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm)   ! linear free surface
                END_2D
             ELSE                                   ! no ocean cavities (only ocean surface)
@@ -356 +358 @@
       !! ** Purpose :  Computation of advective flux with Quickest scheme
       !!
-      !! ** Method :   
+      !! ** Method :
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfu   ! second upwind point
@@ -363 +365 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux
       !!
-      INTEGER  ::  ji, jj, jk               ! dummy loop indices 
-      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars          
+      INTEGER  ::  ji, jj, jk               ! dummy loop indices
+      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars
       REAL(wp) ::  zc, zcurv, zfho          !   -      -
       !----------------------------------------------------------------------
@@ -374 +376 @@
          zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
          zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
-         zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST 
+         zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
          !
          zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
@@ -380 +382 @@
          zcoef3 = ABS( zcurv )
          IF( zcoef3 >= zcoef2 ) THEN
-            zfho = pfc(ji,jj,jk) 
+            zfho = pfc(ji,jj,jk)
          ELSE
             zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
             IF( zcoef1 >= 0. ) THEN
-               zfho = MAX( pfc(ji,jj,jk), zfho ) 
-               zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 
+               zfho = MAX( pfc(ji,jj,jk), zfho )
+               zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) )
             ELSE
-               zfho = MIN( pfc(ji,jj,jk), zfho ) 
-               zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 
+               zfho = MIN( pfc(ji,jj,jk), zfho )
+               zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) )
             ENDIF
          ENDIF

NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traadv_ubs.F90

@@ -10 +10 @@
    !!----------------------------------------------------------------------
    !!   tra_adv_ubs : update the tracer trend with the horizontal
-   !!                 advection trends using a third order biaised scheme  
+   !!                 advection trends using a third order biaised scheme
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and active tracers
@@ -16 +16 @@
    USE trc_oce        ! share passive tracers/Ocean variables
    USE trd_oce        ! trends: ocean variables
-   USE traadv_fct      ! acces to routine interp_4th_cpt 
-   USE trdtra         ! trends manager: tracers 
+   USE traadv_fct      ! acces to routine interp_4th_cpt
+   USE trdtra         ! trends manager: tracers
    USE diaptr         ! poleward transport diagnostics
    USE diaar5         ! AR5 diagnostics
@@ -25 +25 @@
    USE lib_mpp        ! massively parallel library
    USE lbclnk         ! ocean lateral boundary condition (or mpp link)
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
 
    IMPLICIT NONE
@@ -39 +39 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -50 +51 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv_ubs  ***
-      !!                 
+      !!
       !! ** Purpose :   Compute the now trend due to the advection of tracers
       !!      and add it to the general trend of passive tracer equations.
@@ -59 +60 @@
       !!      For example the i-component of the advective fluxes are given by :
       !!                !  e2u e3u un ( mi(Tn) - zltu(i  ) )   if un(i) >= 0
-      !!          ztu = !  or 
+      !!          ztu = !  or
       !!                !  e2u e3u un ( mi(Tn) - zltu(i+1) )   if un(i) < 0
       !!      where zltu is the second derivative of the before temperature field:
       !!          zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ]
-      !!        This results in a dissipatively dominant (i.e. hyper-diffusive) 
-      !!      truncation error. The overall performance of the advection scheme 
-      !!      is similar to that reported in (Farrow and Stevens, 1995). 
+      !!        This results in a dissipatively dominant (i.e. hyper-diffusive)
+      !!      truncation error. The overall performance of the advection scheme
+      !!      is similar to that reported in (Farrow and Stevens, 1995).
       !!        For stability reasons, the first term of the fluxes which corresponds
-      !!      to a second order centered scheme is evaluated using the now velocity 
-      !!      (centered in time) while the second term which is the diffusive part 
-      !!      of the scheme, is evaluated using the before velocity (forward in time). 
+      !!      to a second order centered scheme is evaluated using the now velocity
+      !!      (centered in time) while the second term which is the diffusive part
+      !!      of the scheme, is evaluated using the before velocity (forward in time).
       !!      Note that UBS is not positive. Do not use it on passive tracers.
       !!                On the vertical, the advection is evaluated using a FCT scheme,
-      !!      as the UBS have been found to be too diffusive. 
-      !!                kn_ubs_v argument controles whether the FCT is based on 
-      !!      a 2nd order centrered scheme (kn_ubs_v=2) or on a 4th order compact 
+      !!      as the UBS have been found to be too diffusive.
+      !!                kn_ubs_v argument controles whether the FCT is based on
+      !!      a 2nd order centrered scheme (kn_ubs_v=2) or on a 4th order compact
       !!      scheme (kn_ubs_v=4).
       !!
@@ -81 +82 @@
       !!             - poleward advective heat and salt transport (ln_diaptr=T)
       !!
-      !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404. 
-      !!             Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731Ð1741. 
+      !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404.
+      !!             Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731�1741.
       !!----------------------------------------------------------------------
       INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index
@@ -111 +112 @@
       l_ptr = .FALSE.
       IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
-      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )   l_ptr = .TRUE. 
+      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )   l_ptr = .TRUE.
       IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
          &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE.
@@ -122 +123 @@
       DO jn = 1, kjpt                                            ! tracer loop
          !                                                       ! ===========
-         !                                              
+         !
          DO jk = 1, jpkm1        !==  horizontal laplacian of before tracer ==!
             DO_2D_10_10
@@ -135 +136 @@
                zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef
             END_2D
-            !                                    
-         END DO         
+            !
+         END DO
          CALL lbc_lnk( 'traadv_ubs', zltu, 'T', 1. )   ;    CALL lbc_lnk( 'traadv_ubs', zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn)
-         !    
+         !
          DO_3D_10_10( 1, jpkm1 )
             zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) )      ! upstream transport (x2)
@@ -158 +159 @@
                pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs)                        &
                   &             - (  ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk)    &
-                  &                + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk)  ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+                  &                + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk)  ) &
+                  &                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
             END_2D
-            !                                             
+            !
          END DO
          !
          zltu(:,:,:) = pt(:,:,:,jn,Krhs) - zltu(:,:,:)    ! Horizontal advective trend used in vertical 2nd order FCT case
-         !                                            ! and/or in trend diagnostic (l_trd=T) 
-         !                
+         !                                            ! and/or in trend diagnostic (l_trd=T)
+         !
          IF( l_trd ) THEN                  ! trend diagnostics
              CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, ztu, pU, pt(:,:,:,jn,Kmm) )
              CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, ztv, pV, pt(:,:,:,jn,Kmm) )
          END IF
-         !     
+         !
          !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes)
          IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', ztv(:,:,:) )
@@ -181 +183 @@
          SELECT CASE( kn_ubs_v )       ! select the vertical advection scheme
          !
-         CASE(  2  )                   ! 2nd order FCT 
-            !         
+         CASE(  2  )                   ! 2nd order FCT
+            !
             IF( l_trd )   zltv(:,:,:) = pt(:,:,:,jn,Krhs)          ! store pt(:,:,:,:,Krhs) if trend diag.
             !
@@ -194 +196 @@
                IF( ln_isfcav ) THEN                ! top of the ice-shelf cavities and at the ocean surface
                   DO_2D_11_11
-                     ztw(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface 
+                     ztw(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface
                   END_2D
                ELSE                                ! no cavities: only at the ocean surface
@@ -202 +204 @@
             !
             DO_3D_00_00( 1, jpkm1 )
-               ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
-               pt(ji,jj,jk,jn,Krhs) =   pt(ji,jj,jk,jn,Krhs) +  ztak 
+               ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )    &
+                  &     * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+               pt(ji,jj,jk,jn,Krhs) =   pt(ji,jj,jk,jn,Krhs) +  ztak
                zti(ji,jj,jk)    = ( pt(ji,jj,jk,jn,Kbb) + p2dt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
             END_3D
@@ -228 +231 @@
          !
          DO_3D_00_00( 1, jpkm1 )
-            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )    &
+               &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
          END_3D
          !
@@ -235 +239 @@
                zltv(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs) - zltv(ji,jj,jk)                          &
                   &           + pt(ji,jj,jk,jn,Kmm) * (  pW(ji,jj,jk) - pW(ji,jj,jk+1)  )   &
-                  &                              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+                  &                                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
             END_3D
             CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zltv )
@@ -248 +252 @@
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE nonosc_z  ***
-      !!     
-      !! **  Purpose :   compute monotonic tracer fluxes from the upstream 
-      !!       scheme and the before field by a nonoscillatory algorithm 
+      !!
+      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
+      !!       scheme and the before field by a nonoscillatory algorithm
       !!
       !! **  Method  :   ... ???

NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/trabbc.F90

@@ -12 +12 @@
 
    !!----------------------------------------------------------------------
-   !!   tra_bbc       : update the tracer trend at ocean bottom 
+   !!   tra_bbc       : update the tracer trend at ocean bottom
    !!   tra_bbc_init  : initialization of geothermal heat flux trend
    !!----------------------------------------------------------------------
@@ -19 +19 @@
    USE phycst         ! physical constants
    USE trd_oce        ! trends: ocean variables
-   USE trdtra         ! trends manager: tracers 
+   USE trdtra         ! trends manager: tracers
    !
    USE in_out_manager ! I/O manager
-   USE iom            ! xIOS 
+   USE iom            ! xIOS
    USE fldread        ! read input fields
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
@@ -43 +43 @@
 
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qgh   ! structure of input qgh (file informations, fields read)
- 
+
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -57 +58 @@
       !!                  ***  ROUTINE tra_bbc  ***
       !!
-      !! ** Purpose :   Compute the bottom boundary contition on temperature 
-      !!              associated with geothermal heating and add it to the 
+      !! ** Purpose :   Compute the bottom boundary contition on temperature
+      !!              associated with geothermal heating and add it to the
       !!              general trend of temperature equations.
       !!
-      !! ** Method  :   The geothermal heat flux set to its constant value of 
+      !! ** Method  :   The geothermal heat flux set to its constant value of
       !!              86.4 mW/m2 (Stein and Stein 1992, Huang 1999).
       !!       The temperature trend associated to this heat flux through the
@@ -91 +92 @@
       !                             !  Add the geothermal trend on temperature
       DO_2D_00_00
-         pts(ji,jj,mbkt(ji,jj),jp_tem,Krhs) = pts(ji,jj,mbkt(ji,jj),jp_tem,Krhs) + qgh_trd0(ji,jj) / e3t(ji,jj,mbkt(ji,jj),Kmm)
+         pts(ji,jj,mbkt(ji,jj),jp_tem,Krhs) = pts(ji,jj,mbkt(ji,jj),jp_tem,Krhs)   &
+            &             + qgh_trd0(ji,jj) / e3t(ji,jj,mbkt(ji,jj),Kmm)
       END_2D
       !
@@ -133 +135 @@
       CHARACTER(len=256) ::   cn_dir    ! Root directory for location of ssr files
       !!
-      NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst, sn_qgh, cn_dir 
+      NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst, sn_qgh, cn_dir
       !!----------------------------------------------------------------------
       !

NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/trabbl.F90

@@ -31 +31 @@
    USE trdtra         ! trends: active tracers
    !
-   USE iom            ! IOM library               
+   USE iom            ! IOM library
    USE in_out_manager ! I/O manager
    USE lbclnk         ! ocean lateral boundary conditions
    USE prtctl         ! Print control
    USE timing         ! Timing
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
 
    IMPLICIT NONE
@@ -68 +68 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -195 +196 @@
             zptb(ji,jj) = pt(ji,jj,ik,jn)                ! bottom before T and S
          END_2D
-         !               
+         !
          DO_2D_00_00
             ik = mbkt(ji,jj)                            ! bottom T-level index
@@ -391 +392 @@
                za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
                zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
-               !                                                          ! 2*masked bottom density gradient 
+               !                                                          ! 2*masked bottom density gradient
                zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
                          - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
@@ -513 +514 @@
       END_2D
       ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
-      zmbku(:,:) = REAL( mbku_d(:,:), wp )   ;     zmbkv(:,:) = REAL( mbkv_d(:,:), wp )  
-      CALL lbc_lnk_multi( 'trabbl', zmbku,'U',1., zmbkv,'V',1.) 
+      zmbku(:,:) = REAL( mbku_d(:,:), wp )   ;     zmbkv(:,:) = REAL( mbkv_d(:,:), wp )
+      CALL lbc_lnk_multi( 'trabbl', zmbku,'U',1., zmbkv,'V',1.)
       mbku_d(:,:) = MAX( INT( zmbku(:,:) ), 1 ) ;  mbkv_d(:,:) = MAX( NINT( zmbkv(:,:) ), 1 )
       !

NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traisf.F90

@@ -33 +33 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_isf  ***
-      !!                   
+      !!
       !! ** Purpose :  Compute the temperature trend due to the ice shelf melting (qhoce + qhc)
       !!
@@ -61 +61 @@
          !
          ! Dynamical stability at start up after change in under ice shelf cavity geometry is achieve by correcting the divergence.
-         ! This is achieved by applying a volume flux in order to keep the horizontal divergence after remapping 
+         ! This is achieved by applying a volume flux in order to keep the horizontal divergence after remapping
          ! the same as at the end of the latest time step. So correction need to be apply at nit000 (euler time step) and
          ! half of it at nit000+1 (leap frog time step).
@@ -89 +89 @@
       !! *** Purpose :  Compute the temperature trend due to the ice shelf melting (qhoce + qhc) for cav or par case
       !!
-      !! *** Action :: Update pts(:,:,:,:,Krhs) with the surface boundary condition trend 
+      !! *** Action :: Update pts(:,:,:,:,Krhs) with the surface boundary condition trend
       !!
       !!----------------------------------------------------------------------
@@ -98 +98 @@
       REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) :: ptsc , ptsc_b
       !!----------------------------------------------------------------------
-      INTEGER                      :: ji,jj,jk  ! loop index   
+      INTEGER                      :: ji,jj,jk  ! loop index
       INTEGER                      :: ikt, ikb  ! top and bottom level of the tbl
       REAL(wp), DIMENSION(jpi,jpj) :: ztc       ! total ice shelf tracer trend
@@ -117 +117 @@
          END DO
          !
-         ! level partially include in ice shelf boundary layer 
+         ! level partially include in ice shelf boundary layer
          pts(ji,jj,ikb,jp_tem) = pts(ji,jj,ikb,jp_tem) + ztc(ji,jj) * pfrac(ji,jj)
          !
@@ -128 +128 @@
       !!                  ***  ROUTINE tra_isf_cpl  ***
       !!
-      !! *** Action :: Update pts(:,:,:,:,Krhs) with the ice shelf coupling trend 
+      !! *** Action :: Update pts(:,:,:,:,Krhs) with the ice shelf coupling trend
       !!
       !!----------------------------------------------------------------------

NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traldf_iso.F90

@@ -15 +15 @@
    !!----------------------------------------------------------------------
    !!   tra_ldf_iso   : update the tracer trend with the horizontal component of a iso-neutral laplacian operator
-   !!                   and with the vertical part of the isopycnal or geopotential s-coord. operator 
+   !!                   and with the vertical part of the isopycnal or geopotential s-coord. operator
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and active tracers
@@ -41 +41 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -54 +55 @@
       !!                  ***  ROUTINE tra_ldf_iso  ***
       !!
-      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive 
-      !!      trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and 
+      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
+      !!      trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and
       !!      add it to the general trend of tracer equation.
       !!
-      !! ** Method  :   The horizontal component of the lateral diffusive trends 
+      !! ** Method  :   The horizontal component of the lateral diffusive trends
       !!      is provided by a 2nd order operator rotated along neural or geopo-
       !!      tential surfaces to which an eddy induced advection can be added
@@ -69 +70 @@
       !!
       !!      2nd part :  horizontal fluxes of the lateral mixing operator
-      !!      ========    
+      !!      ========
       !!         zftu =  pahu e2u*e3u/e1u di[ tb ]
       !!               - pahu e2u*uslp    dk[ mi(mk(tb)) ]
@@ -111 +112 @@
       REAL(wp) ::  zcoef0, ze3w_2, zsign, z2dt, z1_2dt   !   -      -
       REAL(wp), DIMENSION(jpi,jpj)     ::   zdkt, zdk1t, z2d
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdit, zdjt, zftu, zftv, ztfw 
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdit, zdjt, zftu, zftv, ztfw
       !!----------------------------------------------------------------------
       !
@@ -119 +120 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
          !
-         akz     (:,:,:) = 0._wp      
+         akz     (:,:,:) = 0._wp
          ah_wslp2(:,:,:) = 0._wp
       ENDIF
-      !   
+      !
       l_hst = .FALSE.
       l_ptr = .FALSE.
-      IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) )     l_ptr = .TRUE. 
+      IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) )     l_ptr = .TRUE.
       IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
          &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE.
@@ -138 +139 @@
       ELSE                    ;   zsign = -1._wp
       ENDIF
-         
+
       !!----------------------------------------------------------------------
       !!   0 - calculate  ah_wslp2 and akz
@@ -173 +174 @@
                DO_3D_10_10( 2, jpkm1 )
                   akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   &
-                     &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) )  )
+                     &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk)    &
+                     &                        / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) )  )
                END_3D
             ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator
@@ -184 +186 @@
            !
          ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit
-            akz(:,:,:) = ah_wslp2(:,:,:)      
+            akz(:,:,:) = ah_wslp2(:,:,:)
          ENDIF
       ENDIF
@@ -191 +193 @@
       DO jn = 1, kjpt                                            ! tracer loop
          !                                                       ! ===========
-         !                                               
-         !!----------------------------------------------------------------------
-         !!   I - masked horizontal derivative 
+         !
+         !!----------------------------------------------------------------------
+         !!   I - masked horizontal derivative
          !!----------------------------------------------------------------------
 !!gm : bug.... why (x,:,:)?   (1,jpj,:) and (jpi,1,:) should be sufficient....
@@ -200 +202 @@
          !!end
 
-         ! Horizontal tracer gradient 
+         ! Horizontal tracer gradient
          DO_3D_10_10( 1, jpkm1 )
             zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)
@@ -207 +209 @@
          IF( ln_zps ) THEN      ! botton and surface ocean correction of the horizontal gradient
             DO_2D_10_10
-               zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)          
+               zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)
                zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
             END_2D
             IF( ln_isfcav ) THEN      ! first wet level beneath a cavity
                DO_2D_10_10
-                  IF( miku(ji,jj) > 1 )   zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)          
-                  IF( mikv(ji,jj) > 1 )   zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)     
+                  IF( miku(ji,jj) > 1 )   zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)
+                  IF( mikv(ji,jj) > 1 )   zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)
                END_2D
             ENDIF
@@ -248 +250 @@
                zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   &
                   &               + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj)      &
-                  &                          + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk)                  
+                  &                          + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk)
             END_2D
             !
@@ -254 +256 @@
                pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      &
                   &                                                 + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   &
-                  &                                              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+                  &                                             * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
             END_2D
-         END DO                                        !   End of slab  
+         END DO                                        !   End of slab
 
          !!----------------------------------------------------------------------
@@ -266 +268 @@
          !                          ! Surface and bottom vertical fluxes set to zero
          ztfw(:,:, 1 ) = 0._wp      ;      ztfw(:,:,jpk) = 0._wp
-         
+
          DO_3D_00_00( 2, jpkm1 )
             !
@@ -295 +297 @@
             END_3D
             !
-         ELSE                                   ! bilaplacian 
+         ELSE                                   ! bilaplacian
             SELECT CASE( kpass )
             CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2
@@ -311 +313 @@
             END SELECT
          ENDIF
-         !         
+         !
          DO_3D_00_00( 1, jpkm1 )
             pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   &

NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traldf_lap_blp.F90

@@ -4 +4 @@
    !! Ocean tracers:  lateral diffusivity trend  (laplacian and bilaplacian)
    !!==============================================================================
-   !! History :  3.7  ! 2014-01  (G. Madec, S. Masson)  Original code, re-entrant laplacian 
+   !! History :  3.7  ! 2014-01  (G. Madec, S. Masson)  Original code, re-entrant laplacian
    !!----------------------------------------------------------------------
 
@@ -38 +38 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -47 +48 @@
    SUBROUTINE tra_ldf_lap( kt, Kmm, kit000, cdtype, pahu, pahv  ,               &
       &                                             pgu , pgv   , pgui, pgvi,   &
-      &                                             pt  , pt_rhs, kjpt, kpass ) 
+      &                                             pt  , pt_rhs, kjpt, kpass )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_ldf_lap  ***
-      !!                   
-      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive 
+      !!
+      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
       !!      trend and add it to the general trend of tracer equation.
       !!
       !! ** Method  :   Second order diffusive operator evaluated using before
-      !!      fields (forward time scheme). The horizontal diffusive trends of 
+      !!      fields (forward time scheme). The horizontal diffusive trends of
       !!      the tracer is given by:
       !!          difft = 1/(e1e2t*e3t) {  di-1[ pahu e2u*e3u/e1u di(tb) ]
@@ -62 +63 @@
       !!          pt_rhs = pt_rhs + difft
       !!
-      !! ** Action  : - Update pt_rhs arrays with the before iso-level 
+      !! ** Action  : - Update pt_rhs arrays with the before iso-level
       !!                harmonic mixing trend.
       !!----------------------------------------------------------------------
@@ -75 +76 @@
       REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt         ! before tracer fields
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs     ! tracer trend 
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs     ! tracer trend
       !
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
@@ -105 +106 @@
       !                             ! =========== !
       DO jn = 1, kjpt               ! tracer loop !
-         !                          ! =========== !    
-         !                               
+         !                          ! =========== !
+         !
          DO_3D_10_10( 1, jpkm1 )
             ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) )
@@ -118 +119 @@
             IF( ln_isfcav ) THEN                ! top in ocean cavities only
                DO_2D_10_10
-                  IF( miku(ji,jj) > 1 )   ztu(ji,jj,miku(ji,jj)) = zaheeu(ji,jj,miku(ji,jj)) * pgui(ji,jj,jn) 
-                  IF( mikv(ji,jj) > 1 )   ztv(ji,jj,mikv(ji,jj)) = zaheev(ji,jj,mikv(ji,jj)) * pgvi(ji,jj,jn) 
+                  IF( miku(ji,jj) > 1 )   ztu(ji,jj,miku(ji,jj)) = zaheeu(ji,jj,miku(ji,jj)) * pgui(ji,jj,jn)
+                  IF( mikv(ji,jj) > 1 )   ztv(ji,jj,mikv(ji,jj)) = zaheev(ji,jj,mikv(ji,jj)) * pgvi(ji,jj,jn)
                END_2D
             ENDIF
@@ -142 +143 @@
       !
    END SUBROUTINE tra_ldf_lap
-   
+
 
    SUBROUTINE tra_ldf_blp( kt, Kmm, kit000, cdtype, pahu, pahv  ,             &
@@ -149 +150 @@
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE tra_ldf_blp  ***
-      !!                    
-      !! ** Purpose :   Compute the before lateral tracer diffusive 
+      !!
+      !! ** Purpose :   Compute the before lateral tracer diffusive
       !!      trend and add it to the general trend of tracer equation.
       !!
@@ -200 +201 @@
       !
       CALL lbc_lnk( 'traldf_lap_blp', zlap(:,:,:,:) , 'T', 1. )     ! Lateral boundary conditions (unchanged sign)
-      !                                               ! Partial top/bottom cell: GRADh( zlap )  
+      !                                               ! Partial top/bottom cell: GRADh( zlap )
       IF( ln_isfcav .AND. ln_zps ) THEN   ;   CALL zps_hde_isf( kt, Kmm, kjpt, zlap, zglu, zglv, zgui, zgvi )  ! both top & bottom
-      ELSEIF(             ln_zps ) THEN   ;   CALL zps_hde    ( kt, Kmm, kjpt, zlap, zglu, zglv )              ! only bottom 
+      ELSEIF(             ln_zps ) THEN   ;   CALL zps_hde    ( kt, Kmm, kjpt, zlap, zglu, zglv )              ! only bottom
       ENDIF
       !

NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traldf_triad.F90

@@ -41 +41 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -108 +109 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
       ENDIF
-      !   
+      !
       l_hst = .FALSE.
       l_ptr = .FALSE.
-      IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) )      l_ptr = .TRUE. 
+      IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) )      l_ptr = .TRUE.
       IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
          &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE.
@@ -124 +125 @@
       ELSE                    ;   zsign = -1._wp
       ENDIF
-      !    
+      !
       !!----------------------------------------------------------------------
       !!   0 - calculate  ah_wslp2, akz, and optionally zpsi_uw, zpsi_vw
@@ -131 +132 @@
       IF( kpass == 1 ) THEN         !==  first pass only  and whatever the tracer is  ==!
          !
-         akz     (:,:,:) = 0._wp      
+         akz     (:,:,:) = 0._wp
          ah_wslp2(:,:,:) = 0._wp
          IF( ln_ldfeiv_dia ) THEN
@@ -158 +159 @@
          END DO
          !
-         DO jp = 0, 1                            ! j-k triads 
+         DO jp = 0, 1                            ! j-k triads
             DO kp = 0, 1
                DO_3D_10_10( 1, jpkm1 )
@@ -184 +185 @@
                DO_3D_10_10( 2, jpkm1 )
                   akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   &
-                     &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) )  )
+                     &          * (  akz(ji,jj,jk)              &
+                     &            + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) )  )
                END_3D
             ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator
@@ -195 +197 @@
            !
          ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit
-            akz(:,:,:) = ah_wslp2(:,:,:)      
+            akz(:,:,:) = ah_wslp2(:,:,:)
          ENDIF
          !
@@ -222 +224 @@
             IF( ln_isfcav ) THEN                   ! top level (ocean cavities only)
                DO_2D_10_10
-                  IF( miku(ji,jj)  > 1 )   zdit(ji,jj,miku(ji,jj) ) = pgui(ji,jj,jn) 
-                  IF( mikv(ji,jj)  > 1 )   zdjt(ji,jj,mikv(ji,jj) ) = pgvi(ji,jj,jn) 
+                  IF( miku(ji,jj)  > 1 )   zdit(ji,jj,miku(ji,jj) ) = pgui(ji,jj,jn)
+                  IF( mikv(ji,jj)  > 1 )   zdjt(ji,jj,mikv(ji,jj) ) = pgvi(ji,jj,jn)
                END_2D
             ENDIF
@@ -330 +332 @@
             !                             !==  horizontal divergence and add to the general trend  ==!
             DO_2D_00_00
-               pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  zftu(ji-1,jj,jk) - zftu(ji,jj,jk)       &
-                  &                                           + zftv(ji,jj-1,jk) - zftv(ji,jj,jk)   )   &
-                  &                                        / (  e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm)  )
+               pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    &
+                  &                            + zsign * (  zftu(ji-1,jj,jk) - zftu(ji,jj,jk)       &
+                  &                                       + zftv(ji,jj-1,jk) - zftv(ji,jj,jk)   )   &
+                  &                                     / (  e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm)  )
             END_2D
             !
@@ -344 +347 @@
                   &                            * (  pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )
             END_3D
-         ELSE                                   ! bilaplacian 
+         ELSE                                   ! bilaplacian
             SELECT CASE( kpass )
             CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2
@@ -357 +360 @@
                      &                               + akz     (ji,jj,jk) * ( pt2(ji,jj,jk-1,jn) - pt2(ji,jj,jk,jn) )   )
                END_3D
-            END SELECT 
+            END SELECT
          ENDIF
          !
          DO_3D_00_00( 1, jpkm1 )
-            pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   &
-               &                                              / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) )
+            pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    &
+            &                                  + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   &
+               &                                       / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) )
          END_3D
          !

NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/tramle.F90

@@ -49 +49 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)

NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/tranpc.F90

@@ -35 +35 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -71 +72 @@
       REAL(wp), DIMENSION(        jpk     )   ::   zvn2         ! vertical profile of N2 at 1 given point...
       REAL(wp), DIMENSION(        jpk,jpts)   ::   zvts, zvab   ! vertical profile of T & S , and  alpha & betaat 1 given point
-      REAL(wp), DIMENSION(jpi,jpj,jpk     )   ::   zn2          ! N^2 
+      REAL(wp), DIMENSION(jpi,jpj,jpk     )   ::   zn2          ! N^2
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts)   ::   zab          ! alpha and beta
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds ! 3D workspace
@@ -86 +87 @@
          IF( l_trdtra )   THEN                    !* Save initial after fields
             ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
-            ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kaa) 
+            ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kaa)
             ztrds(:,:,:) = pts(:,:,:,jp_sal,Kaa)
          ENDIF
@@ -92 +93 @@
          IF( l_LB_debug ) THEN
             ! Location of 1 known convection site to follow what's happening in the water column
-            ilc1 = 45 ;  jlc1 = 3 ; !  ORCA2 4x4, Antarctic coast, more than 2 unstable portions in the  water column...           
+            ilc1 = 45 ;  jlc1 = 3 ; !  ORCA2 4x4, Antarctic coast, more than 2 unstable portions in the  water column...
             nncpu = 1  ;            ! the CPU domain contains the convection spot
-            klc1 =  mbkt(ilc1,jlc1) ! bottom of the ocean for debug point...          
+            klc1 =  mbkt(ilc1,jlc1) ! bottom of the ocean for debug point...
          ENDIF
          !
@@ -105 +106 @@
             !
             IF( tmask(ji,jj,2) == 1 ) THEN      ! At least 2 ocean points
-               !                                     ! consider one ocean column 
+               !                                     ! consider one ocean column
                zvts(:,jp_tem) = pts(ji,jj,:,jp_tem,Kaa)      ! temperature
                zvts(:,jp_sal) = pts(ji,jj,:,jp_sal,Kaa)      ! salinity
                !
-               zvab(:,jp_tem)  = zab(ji,jj,:,jp_tem)     ! Alpha 
-               zvab(:,jp_sal)  = zab(ji,jj,:,jp_sal)     ! Beta  
-               zvn2(:)         = zn2(ji,jj,:)            ! N^2 
+               zvab(:,jp_tem)  = zab(ji,jj,:,jp_tem)     ! Alpha
+               zvab(:,jp_sal)  = zab(ji,jj,:,jp_sal)     ! Beta
+               zvn2(:)         = zn2(ji,jj,:)            ! N^2
                !
                IF( l_LB_debug ) THEN                  !LB debug:
@@ -117 +118 @@
                   IF( ( ji == ilc1 ).AND.( jj == jlc1 ) ) lp_monitor_point = .TRUE.
                   ! writing only if on CPU domain where conv region is:
-                  lp_monitor_point = (narea == nncpu).AND.lp_monitor_point                      
+                  lp_monitor_point = (narea == nncpu).AND.lp_monitor_point
                ENDIF                                  !LB debug  end
                !
@@ -129 +130 @@
                   !
                   jiter = jiter + 1
-                  ! 
+                  !
                   IF( jiter >= 400 ) EXIT
                   !
@@ -144 +145 @@
                         ilayer = ilayer + 1    ! yet another instable portion of the water column found....
                         !
-                        IF( lp_monitor_point ) THEN 
+                        IF( lp_monitor_point ) THEN
                            WRITE(numout,*)
                            IF( ilayer == 1 .AND. jiter == 1 ) THEN   ! first time a column is spoted with an instability
@@ -159 +160 @@
                         ENDIF
                         !
-                        IF( jiter == 1 )   inpcc = inpcc + 1 
+                        IF( jiter == 1 )   inpcc = inpcc + 1
                         !
                         IF( lp_monitor_point )   WRITE(numout, *) 'Negative N2 at ikp =',ikp,' for layer #', ilayer
@@ -184 +185 @@
                         zsum_beta = 0._wp
                         zsum_z    = 0._wp
-                                                 
+
                         DO jk = ikup, ikbot      ! Inside the instable (and overlying neutral) portion of the column
                            !
@@ -193 +194 @@
                            zsum_beta = zsum_beta + zvab(jk,jp_sal)*zdz
                            zsum_z    = zsum_z    + zdz
-                           !                              
+                           !
                            IF( jk == ikbot ) EXIT ! avoid array-index overshoot in case ikbot = jpk, cause we're calling jk+1 next line
                            !! EXIT when we have reached the last layer that is instable (N2<0) or neutral (N2=0):
                            IF( zvn2(jk+1) > zn2_zero ) EXIT
                         END DO
-                       
+
                         ikdown = jk ! for the current unstable layer, ikdown is the deepest point with a negative or neutral N2
                         IF( ikup == ikdown )   CALL ctl_stop( 'tra_npc :  PROBLEM #2')
@@ -224 +225 @@
                            zvab(jk,jp_sal) = zbeta
                         END DO
-                        
-                        
+
+
                         !! Updating N2 in the relvant portion of the water column
                         !! Temperature, Salinity, Alpha and Beta have been homogenized in the unstable portion
                         !! => Need to re-compute N2! will use Alpha and Beta!
-                        
+
                         ikup   = MAX(2,ikup)         ! ikup can never be 1 !
                         ik_low = MIN(ikdown+1,ikbot) ! we must go 1 point deeper than ikdown!
-                        
+
                         DO jk = ikup, ik_low              ! we must go 1 point deeper than ikdown!
 
@@ -252 +253 @@
 
                         END DO
-                     
+
                         ikp = MIN(ikdown+1,ikbot)
-                        
+
 
                      ENDIF  !IF( zvn2(ikp) < 0. )
@@ -264 +265 @@
 
                   IF( ikp /= ikbot )   CALL ctl_stop( 'tra_npc :  PROBLEM #3')
-                 
+
                   ! ******* At this stage ikp == ikbot ! *******
-                 
+
                   IF( ilayer > 0 ) THEN      !! least an unstable layer has been found
                      !

NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traqsr.F90

@@ -9 +9 @@
    !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
    !!             -   !  2005-11  (G. Madec) zco, zps, sco coordinate
-   !!            3.2  !  2009-04  (G. Madec & NEMO team) 
-   !!            3.6  !  2012-05  (C. Rousset) store attenuation coef for use in ice model 
+   !!            3.2  !  2009-04  (G. Madec & NEMO team)
+   !!            3.6  !  2012-05  (C. Rousset) store attenuation coef for use in ice model
    !!            3.6  !  2015-12  (O. Aumont, J. Jouanno, C. Ethe) use vertical profile of chlorophyll
-   !!            3.7  !  2015-11  (G. Madec, A. Coward)  remove optimisation for fix volume 
+   !!            3.7  !  2015-11  (G. Madec, A. Coward)  remove optimisation for fix volume
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
-   !!   tra_qsr       : temperature trend due to the penetration of solar radiation 
-   !!   tra_qsr_init  : initialization of the qsr penetration 
+   !!   tra_qsr       : temperature trend due to the penetration of solar radiation
+   !!   tra_qsr_init  : initialization of the qsr penetration
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and active tracers
@@ -44 +44 @@
    !                                 !!* Namelist namtra_qsr: penetrative solar radiation
    LOGICAL , PUBLIC ::   ln_traqsr    !: light absorption (qsr) flag
-   LOGICAL , PUBLIC ::   ln_qsr_rgb   !: Red-Green-Blue light absorption flag  
+   LOGICAL , PUBLIC ::   ln_qsr_rgb   !: Red-Green-Blue light absorption flag
    LOGICAL , PUBLIC ::   ln_qsr_2bd   !: 2 band         light absorption flag
    LOGICAL , PUBLIC ::   ln_qsr_bio   !: bio-model      light absorption flag
@@ -53 +53 @@
    !
    INTEGER , PUBLIC ::   nksr         !: levels below which the light cannot penetrate (depth larger than 391 m)
- 
+
    INTEGER, PARAMETER ::   np_RGB  = 1   ! R-G-B     light penetration with constant Chlorophyll
    INTEGER, PARAMETER ::   np_RGBc = 2   ! R-G-B     light penetration with Chlorophyll data
@@ -68 +68 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -86 +87 @@
       !!      Considering the 2 wavebands case:
       !!         I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) )
-      !!         The temperature trend associated with the solar radiation penetration 
+      !!         The temperature trend associated with the solar radiation penetration
       !!         is given by : zta = 1/e3t dk[ I ] / (rau0*Cp)
       !!         At the bottom, boudary condition for the radiation is no flux :
       !!      all heat which has not been absorbed in the above levels is put
       !!      in the last ocean level.
-      !!         The computation is only done down to the level where 
-      !!      I(k) < 1.e-15 W/m2 (i.e. over the top nksr levels) . 
+      !!         The computation is only done down to the level where
+      !!      I(k) < 1.e-15 W/m2 (i.e. over the top nksr levels) .
       !!
       !! ** Action  : - update ta with the penetrative solar radiation trend
@@ -112 +113 @@
       REAL(wp) ::   zz0 , zz1                !    -         -
       REAL(wp) ::   zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze
-      REAL(wp) ::   zlogc, zlogc2, zlogc3 
+      REAL(wp) ::   zlogc, zlogc2, zlogc3
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   :: zekb, zekg, zekr
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt
@@ -127 +128 @@
       !
       IF( l_trdtra ) THEN      ! trends diagnostic: save the input temperature trend
-         ALLOCATE( ztrdt(jpi,jpj,jpk) ) 
+         ALLOCATE( ztrdt(jpi,jpj,jpk) )
          ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
       ENDIF
@@ -163 +164 @@
          ALLOCATE( zekb(jpi,jpj)     , zekg(jpi,jpj)     , zekr  (jpi,jpj)     , &
             &      ze0 (jpi,jpj,jpk) , ze1 (jpi,jpj,jpk) , ze2   (jpi,jpj,jpk) , &
-            &      ze3 (jpi,jpj,jpk) , zea (jpi,jpj,jpk) , zchl3d(jpi,jpj,jpk)   ) 
+            &      ze3 (jpi,jpj,jpk) , zea (jpi,jpj,jpk) , zchl3d(jpi,jpj,jpk)   )
          !
          IF( nqsr == np_RGBc ) THEN          !*  Variable Chlorophyll
@@ -183 +184 @@
                      zpsimax = 0.6   - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3
                      zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2
-                     zCze    = 1.12  * (zchl)**0.803 
+                     zCze    = 1.12  * (zchl)**0.803
                      !
                      zchl3d(ji,jj,jk) = zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) / zdelpsi )**2 ) )
@@ -192 +193 @@
          ELSE                                !* constant chrlorophyll
            DO jk = 1, nksr + 1
-              zchl3d(:,:,jk) = 0.05 
+              zchl3d(:,:,jk) = 0.05
             ENDDO
          ENDIF
@@ -231 +232 @@
          END_3D
          !
-         DEALLOCATE( zekb , zekg , zekr , ze0 , ze1 , ze2 , ze3 , zea , zchl3d ) 
+         DEALLOCATE( zekb , zekg , zekr , ze0 , ze1 , ze2 , ze3 , zea , zchl3d )
          !
       CASE( np_2BD  )            !==  2-bands fluxes  ==!
@@ -240 +241 @@
             zc0 = zz0 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi1r )
             zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi1r )
-            qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) ) 
+            qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) )
          END_3D
          !
@@ -248 +249 @@
       DO_3D_00_00( 1, nksr )
          pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
-            &                      + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t(ji,jj,jk,Kmm)
+            &                      + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) )   &
+            &                             / e3t(ji,jj,jk,Kmm)
       END_3D
       !
@@ -264 +266 @@
          DO jk = nksr, 1, -1
             zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * rau0_rcp
-         END DO         
+         END DO
          CALL iom_put( 'qsr3d', zetot )   ! 3D distribution of shortwave Radiation
-         DEALLOCATE( zetot ) 
+         DEALLOCATE( zetot )
       ENDIF
       !
@@ -272 +274 @@
          IF( lwxios ) CALL iom_swap(      cwxios_context          )
          CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc     , ldxios = lwxios )
-         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev, ldxios = lwxios ) 
+         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev, ldxios = lwxios )
          IF( lwxios ) CALL iom_swap(      cxios_context          )
       ENDIF
@@ -279 +281 @@
          ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt )
-         DEALLOCATE( ztrdt ) 
+         DEALLOCATE( ztrdt )
       ENDIF
       !                       ! print mean trends (used for debugging)
@@ -298 +300 @@
       !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio
       !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The
-      !!      default values correspond to clear water (type I in Jerlov' 
+      !!      default values correspond to clear water (type I in Jerlov'
       !!      (1968) classification.
       !!         called by tra_qsr at the first timestep (nit000)
@@ -348 +350 @@
          &                               ' 2 bands, 3 RGB bands or bio-model light penetration' )
       !
-      IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr = np_RGB 
+      IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr = np_RGB
       IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr = np_RGBc
       IF( ln_qsr_2bd                      )   nqsr = np_2BD
@@ -358 +360 @@
       !
       SELECT CASE( nqsr )
-      !                               
+      !
       CASE( np_RGB , np_RGBc )         !==  Red-Green-Blue light penetration  ==!
-         !                             
+         !
          IF(lwp)   WRITE(numout,*) '   ==>>>   R-G-B   light penetration '
          !
          CALL trc_oce_rgb( rkrgb )                 ! tabulated attenuation coef.
-         !                                   
+         !
          nksr = trc_oce_ext_lev( r_si2, 33._wp )   ! level of light extinction
          !

NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/trasbc.F90

@@ -9 +9 @@
    !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
    !!             -   !  2010-09  (C. Ethe, G. Madec) Merge TRA-TRC
-   !!            3.6  !  2014-11  (P. Mathiot) isf melting forcing 
+   !!            3.6  !  2014-11  (P. Mathiot) isf melting forcing
    !!            4.1  !  2019-09  (P. Mathiot) isf moved in traisf
    !!----------------------------------------------------------------------
@@ -21 +21 @@
    USE phycst         ! physical constant
    USE eosbn2         ! Equation Of State
-   USE sbcmod         ! ln_rnf  
-   USE sbcrnf         ! River runoff  
+   USE sbcmod         ! ln_rnf
+   USE sbcrnf         ! River runoff
    USE traqsr         ! solar radiation penetration
    USE trd_oce        ! trends: ocean variables
-   USE trdtra         ! trends manager: tracers 
-#if defined key_asminc   
+   USE trdtra         ! trends manager: tracers
+#if defined key_asminc
    USE asminc         ! Assimilation increment
 #endif
@@ -43 +43 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -53 +54 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_sbc  ***
-      !!                   
+      !!
       !! ** Purpose :   Compute the tracer surface boundary condition trend of
       !!      (flux through the interface, concentration/dilution effect)
       !!      and add it to the general trend of tracer equations.
       !!
-      !! ** Method :   The (air+ice)-sea flux has two components: 
-      !!      (1) Fext, external forcing (i.e. flux through the (air+ice)-sea interface); 
-      !!      (2) Fwe , tracer carried with the water that is exchanged with air+ice. 
+      !! ** Method :   The (air+ice)-sea flux has two components:
+      !!      (1) Fext, external forcing (i.e. flux through the (air+ice)-sea interface);
+      !!      (2) Fwe , tracer carried with the water that is exchanged with air+ice.
       !!               The input forcing fields (emp, rnf, sfx) contain Fext+Fwe,
       !!             they are simply added to the tracer trend (ts(Krhs)).
@@ -68 +69 @@
       !!             concentration/dilution effect associated with water exchanges.
       !!
-      !! ** Action  : - Update ts(Krhs) with the surface boundary condition trend 
+      !! ** Action  : - Update ts(Krhs) with the surface boundary condition trend
       !!              - send trends to trdtra module for further diagnostics(l_trdtra=T)
       !!----------------------------------------------------------------------
@@ -75 +76 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts        ! active tracers and RHS of tracer equation
       !
-      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices  
+      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
       INTEGER  ::   ikt, ikb                    ! local integers
       REAL(wp) ::   zfact, z1_e3t, zdep, ztim   ! local scalar
@@ -90 +91 @@
       !
       IF( l_trdtra ) THEN                    !* Save ta and sa trends
-         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
+         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
          ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
          ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs)
@@ -127 +128 @@
          sbc_tsc(ji,jj,jp_sal) = r1_rau0     * sfx(ji,jj)   ! salt flux due to freezing/melting
       END_2D
-      IF( ln_linssh ) THEN                !* linear free surface  
+      IF( ln_linssh ) THEN                !* linear free surface
          DO_2D_01_00
             sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rau0 * emp(ji,jj) * pts(ji,jj,1,jp_tem,Kmm)
@@ -138 +139 @@
       DO jn = 1, jpts               !==  update tracer trend  ==!
          DO_2D_01_00
-            pts(ji,jj,1,jn,Krhs) = pts(ji,jj,1,jn,Krhs) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) / e3t(ji,jj,1,Kmm)
+            pts(ji,jj,1,jn,Krhs) = pts(ji,jj,1,jn,Krhs) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) )    &
+               &                                                / e3t(ji,jj,1,Kmm)
          END_2D
       END DO
-      !                  
+      !
       IF( lrst_oce ) THEN           !==  write sbc_tsc in the ocean restart file  ==!
          IF( lwxios ) CALL iom_swap(      cwxios_context          )
@@ -153 +155 @@
       !----------------------------------------
       !
-      IF( ln_rnf ) THEN         ! input of heat and salt due to river runoff 
+      IF( ln_rnf ) THEN         ! input of heat and salt due to river runoff
          zfact = 0.5_wp
          DO_2D_01_00
@@ -162 +164 @@
                                         &                      +  ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep
                   IF( ln_rnf_sal )   pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)                                  &
-                                        &                      +  ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep 
+                                        &                      +  ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep
                END DO
             ENDIF
@@ -179 +181 @@
       IF( ln_sshinc ) THEN         ! input of heat and salt due to assimilation
           !
-         IF( ln_linssh ) THEN 
+         IF( ln_linssh ) THEN
             DO_2D_01_00
                ztim = ssh_iau(ji,jj) / e3t(ji,jj,1,Kmm)
@@ -202 +204 @@
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_nsr, ztrdt )
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_nsr, ztrds )
-         DEALLOCATE( ztrdt , ztrds ) 
+         DEALLOCATE( ztrdt , ztrds )
       ENDIF
       !

NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/trazdf.F90

@@ -13 +13 @@
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers variables
-   USE dom_oce        ! ocean space and time domain variables 
+   USE dom_oce        ! ocean space and time domain variables
    USE domvvl         ! variable volume
    USE phycst         ! physical constant
@@ -19 +19 @@
    USE sbc_oce        ! surface boundary condition: ocean
    USE ldftra         ! lateral diffusion: eddy diffusivity
-   USE ldfslp         ! lateral diffusion: iso-neutral slope 
+   USE ldfslp         ! lateral diffusion: iso-neutral slope
    USE trd_oce        ! trends: ocean variables
    USE trdtra         ! trends: tracer trend manager
@@ -37 +37 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -77 +78 @@
       !
       !                                      !* compute lateral mixing trend and add it to the general trend
-      CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, Kbb, Kmm, Krhs, pts, Kaa, jpts ) 
+      CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, Kbb, Kmm, Krhs, pts, Kaa, jpts )
 
 !!gm WHY here !   and I don't like that !
@@ -88 +89 @@
       IF( l_trdtra )   THEN                      ! save the vertical diffusive trends for further diagnostics
          DO jk = 1, jpkm1
-            ztrdt(:,:,jk) = ( ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_tem,Kbb)*e3t(:,:,jk,Kbb) ) &
-               &          / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrdt(:,:,jk)
-            ztrds(:,:,jk) = ( ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_sal,Kbb)*e3t(:,:,jk,Kbb) ) &
-              &           / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrds(:,:,jk)
+            ztrdt(:,:,jk) = ( ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa)    &
+               &              - pts(:,:,jk,jp_tem,Kbb)*e3t(:,:,jk,Kbb) )  &
+               &             / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrdt(:,:,jk)
+            ztrds(:,:,jk) = ( ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa)    &
+               &              - pts(:,:,jk,jp_sal,Kbb)*e3t(:,:,jk,Kbb) ) &
+               &             / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrds(:,:,jk)
          END DO
 !!gm this should be moved in trdtra.F90 and done on all trends
@@ -108 +111 @@
    END SUBROUTINE tra_zdf
 
- 
-   SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, Kbb, Kmm, Krhs, pt, Kaa, kjpt ) 
+
+   SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, Kbb, Kmm, Krhs, pt, Kaa, kjpt )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_zdf_imp  ***
       !!
       !! ** Purpose :   Compute the after tracer through a implicit computation
-      !!     of the vertical tracer diffusion (including the vertical component 
-      !!     of lateral mixing (only for 2nd order operator, for fourth order 
-      !!     it is already computed and add to the general trend in traldf) 
+      !!     of the vertical tracer diffusion (including the vertical component
+      !!     of lateral mixing (only for 2nd order operator, for fourth order
+      !!     it is already computed and add to the general trend in traldf)
       !!
       !! ** Method  :  The vertical diffusion of a tracer ,t , is given by:
@@ -158 +161 @@
             zwt(:,:,1) = 0._wp
             !
-            IF( l_ldfslp ) THEN            ! isoneutral diffusion: add the contribution 
-               IF( ln_traldf_msc  ) THEN     ! MSC iso-neutral operator 
+            IF( l_ldfslp ) THEN            ! isoneutral diffusion: add the contribution
+               IF( ln_traldf_msc  ) THEN     ! MSC iso-neutral operator
                   DO_3D_00_00( 2, jpkm1 )
-                     zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk)  
+                     zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk)
                   END_3D
                ELSE                          ! standard or triad iso-neutral operator
@@ -204 +207 @@
             !   The solution will be in the 4d array pta.
             !   The 3d array zwt is used as a work space array.
-            !   En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then 
+            !   En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then
             !   used as a work space array: its value is modified.
             !
@@ -214 +217 @@
             END_3D
             !
-         ENDIF 
-         !         
+         ENDIF
+         !
          DO_2D_00_00
-            pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs)
+            pt(ji,jj,1,jn,Kaa) =        e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb)    &
+               &               + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs)
          END_2D
          DO_3D_00_00( 2, jpkm1 )
-            zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs)   ! zrhs=right hand side
+            zrhs =        e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb)    & 
+               & + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs)   ! zrhs=right hand side
             pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa)
          END_3D