New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 5758 for branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Ignore:
Timestamp:
2015-09-24T08:31:40+02:00 (9 years ago)
Author:
gm
Message:

#1593: LDF-ADV, step II.1: phasing the improvements/simplifications of diffusive trend (see wiki)

Location:
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA
Files:
1 added
3 deleted
8 edited
1 moved

Legend:

Unmodified
Added
Removed

  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90

    r5147 r5758  
    66   !! History :  2.0  !  2005-11  (G. Madec)  Original code 
    77   !!            3.3  !  2010-09  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport 
    8    !!            4.0  !  2011-06  (G. Madec)  Addition of Mixed Layer Eddy parameterisation 
     8   !!            3.6  !  2011-06  (G. Madec)  Addition of Mixed Layer Eddy parameterisation 
     9   !!            4.0  !  2014-05  (G. Madec)  Add 2nd/4th order cases for CEN and FCT schemes  
     10   !!             -   !  2014-12  (G. Madec) suppression of cross land advection option 
    911   !!---------------------------------------------------------------------- 
    1012 
    … …  
    2224   USE traadv_ubs      ! UBS      scheme           (tra_adv_ubs    routine) 
    2325   USE traadv_qck      ! QUICKEST scheme           (tra_adv_qck    routine) 
    24    USE traadv_eiv      ! eddy induced velocity     (tra_adv_eiv    routine) 
    2526   USE traadv_mle      ! ML eddy induced velocity  (tra_adv_mle    routine) 
    2627   USE cla             ! cross land advection      (cla_traadv     routine) 
    27    USE ldftra_oce      ! lateral diffusion coefficient on tracers 
     28   USE ldftra          ! lateral diffusion coefficient on tracers 
     29   USE ldfslp          ! Lateral diffusion: slopes of neutral surfaces 
    2830   ! 
    2931   USE in_out_manager  ! I/O manager 
    … …  
    7476      !! ** Method  : - Update (ua,va) with the advection term following nadv 
    7577      !!---------------------------------------------------------------------- 
    76       ! 
    7778      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    7879      ! 
    … …  
    8485      ! 
    8586      CALL wrk_alloc( jpi, jpj, jpk, zun, zvn, zwn ) 
     87      ! 
    8688      !                                          ! set time step 
    8789      IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000 
    … …  
    9597      !                                               !==  effective transport  ==! 
    9698      DO jk = 1, jpkm1 
    97          zun(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)                  ! eulerian transport only 
    98          zvn(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 
    99          zwn(:,:,jk) = e1t(:,:) * e2t(:,:)      * wn(:,:,jk) 
     99         zun(:,:,jk) = e2u  (:,:) * fse3u(:,:,jk) * un(:,:,jk)                  ! eulerian transport only 
     100         zvn(:,:,jk) = e1v  (:,:) * fse3v(:,:,jk) * vn(:,:,jk) 
     101         zwn(:,:,jk) = e1e2t(:,:)                 * wn(:,:,jk) 
    100102      END DO 
    101103      ! 
    … …  
    109111      zwn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom 
    110112      ! 
    111       IF( lk_traldf_eiv .AND. .NOT. ln_traldf_grif )   & 
    112          &              CALL tra_adv_eiv( kt, nit000, zun, zvn, zwn, 'TRA' )    ! add the eiv transport (if necessary) 
    113       ! 
    114       IF( ln_mle    )   CALL tra_adv_mle( kt, nit000, zun, zvn, zwn, 'TRA' )    ! add the mle transport (if necessary) 
    115       ! 
    116       CALL iom_put( "uocetr_eff", zun )                                         ! output effective transport       
     113      IF( ln_ldfeiv .AND. .NOT. ln_traldf_triad )   & 
     114         &              CALL ldf_eiv_trp( kt, nit000, zun, zvn, zwn, 'TRA' )   ! add the eiv transport (if necessary) 
     115      ! 
     116      IF( ln_mle    )   CALL tra_adv_mle( kt, nit000, zun, zvn, zwn, 'TRA' )   ! add the mle transport (if necessary) 
     117      ! 
     118      CALL iom_put( "uocetr_eff", zun )                                        ! output effective transport       
    117119      CALL iom_put( "vocetr_eff", zvn ) 
    118120      CALL iom_put( "wocetr_eff", zwn ) 
    119121      ! 
    120       IF( ln_diaptr )   CALL dia_ptr( zvn )                                     ! diagnose the effective MSF  
     122      IF( ln_diaptr )   CALL dia_ptr( zvn )                                    ! diagnose the effective MSF  
    121123      ! 
    122124    

  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90

    r5120 r5758  
    44   !! Ocean Active tracers : lateral diffusive trends  
    55   !!===================================================================== 
    6    !! History :  9.0  ! 2005-11 (G. Madec)  Original code 
    7    !!       NEMO 3.0  ! 2008-01  (C. Ethe, G. Madec)  merge TRC-TRA  
     6   !! History :  9.0  ! 2005-11  (G. Madec)  Original code 
     7   !!  NEMO      3.0  ! 2008-01  (C. Ethe, G. Madec)  merge TRC-TRA  
     8   !!            3.7  ! 2013-12  (G. Madec) remove the optional computation from T & S anomaly profiles and traldf_bilapg 
     9   !!             -   ! 2013-12  (F. Lemarie, G. Madec)  triad operator (Griffies) + Method of Stabilizing Correction 
     10   !!             -   ! 2014-01  (G. Madec, S. Masson)  restructuration/simplification of lateral diffusive operators 
    811   !!---------------------------------------------------------------------- 
    912 
    … …  
    1114   !!   tra_ldf      : update the tracer trend with the lateral diffusion 
    1215   !!   tra_ldf_init : initialization, namelist read, and parameters control 
    13    !!       ldf_ano  : compute lateral diffusion for constant T-S profiles 
    14    !!---------------------------------------------------------------------- 
    15    USE oce             ! ocean dynamics and tracers 
    16    USE dom_oce         ! ocean space and time domain 
    17    USE phycst          ! physical constants 
    18    USE ldftra_oce      ! ocean tracer   lateral physics 
    19    USE ldfslp          ! ??? 
    20    USE traldf_bilapg   ! lateral mixing            (tra_ldf_bilapg routine) 
    21    USE traldf_bilap    ! lateral mixing             (tra_ldf_bilap routine) 
    22    USE traldf_iso      ! lateral mixing               (tra_ldf_iso routine) 
    23    USE traldf_iso_grif ! lateral mixing          (tra_ldf_iso_grif routine) 
    24    USE traldf_lap      ! lateral mixing               (tra_ldf_lap routine) 
    25    USE trd_oce         ! trends: ocean variables 
    26    USE trdtra          ! trends manager: tracers  
     16   !!---------------------------------------------------------------------- 
     17   USE oce           ! ocean dynamics and tracers 
     18   USE dom_oce       ! ocean space and time domain 
     19   USE phycst        ! physical constants 
     20   USE ldftra        ! lateral diffusion: eddy diffusivity & EIV coeff. 
     21   USE ldfslp        ! lateral diffusion: iso-neutral slope 
     22   USE traldf_lap    ! lateral diffusion: laplacian iso-level            operator  (tra_ldf_lap   routine) 
     23   USE traldf_iso    ! lateral diffusion: laplacian iso-neutral standard operator  (tra_ldf_iso   routine) 
     24   USE traldf_triad  ! lateral diffusion: laplacian iso-neutral triad    operator  (tra_ldf_triad routine) 
     25   USE traldf_blp    ! lateral diffusion (iso-level lap/blp)                       (tra_ldf_lap   routine) 
     26   USE trd_oce       ! trends: ocean variables 
     27   USE trdtra        ! ocean active tracers trends 
    2728   ! 
    28    USE prtctl          ! Print control 
    29    USE in_out_manager  ! I/O manager 
    30    USE lib_mpp         ! distribued memory computing library 
    31    USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    32    USE wrk_nemo        ! Memory allocation 
    33    USE timing          ! Timing 
     29   USE prtctl         ! Print control 
     30   USE in_out_manager ! I/O manager 
     31   USE lib_mpp        ! distribued memory computing library 
     32   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     33   USE wrk_nemo       ! Memory allocation 
     34   USE timing         ! Timing 
    3435 
    3536   IMPLICIT NONE 
    … …  
    3738 
    3839   PUBLIC   tra_ldf        ! called by step.F90  
    39    PUBLIC   tra_ldf_init   ! called by opa.F90  
     40   PUBLIC   tra_ldf_init   ! called by nemogcm.F90  
    4041   ! 
    4142   INTEGER ::   nldf = 0   ! type of lateral diffusion used defined from ln_traldf_... namlist logicals) 
    42  
    43    REAL, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   t0_ldf, s0_ldf   !: lateral diffusion trends of T & S for a cst profile 
    44    !                                                               !  (key_traldf_ano only) 
    45  
     43    
    4644   !! * Substitutions 
    4745#  include "domzgr_substitute.h90" 
    4846#  include "vectopt_loop_substitute.h90" 
    4947   !!---------------------------------------------------------------------- 
    50    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     48   !! NEMO/OPA 3.7 , NEMO Consortium (2015) 
    5149   !! $Id$  
    5250   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    6563      !!---------------------------------------------------------------------- 
    6664      ! 
    67       IF( nn_timing == 1 )  CALL timing_start('tra_ldf') 
    68       ! 
    69       rldf = 1     ! For active tracers the  
    70  
     65      IF( nn_timing == 1 )   CALL timing_start('tra_ldf') 
     66      ! 
    7167      IF( l_trdtra )   THEN                    !* Save ta and sa trends 
    72          CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )  
     68         CALL wrk_alloc( jpi,jpj,jpk,  ztrdt, ztrds )  
    7369         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)  
    7470         ztrds(:,:,:) = tsa(:,:,:,jp_sal) 
    7571      ENDIF 
    76  
    77       SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend 
    78       CASE ( 0 )   ;   CALL tra_ldf_lap     ( kt, nit000, 'TRA', gtsu, gtsv, gtui, gtvi,        & 
    79                                &                                   tsb, tsa, jpts        )  ! iso-level laplacian 
    80       CASE ( 1 )                                                                              ! rotated laplacian 
    81          IF( ln_traldf_grif ) THEN                                                           
    82                        CALL tra_ldf_iso_grif( kt, nit000,'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 )      ! Griffies operator 
    83          ELSE                                                                                 
    84                        CALL tra_ldf_iso     ( kt, nit000, 'TRA', gtsu, gtsv, gtui, gtvi,        & 
    85                                &                                  tsb, tsa, jpts, ahtb0 )      ! Madec operator 
    86          ENDIF 
    87       CASE ( 2 )   ;   CALL tra_ldf_bilap   ( kt, nit000, 'TRA', gtsu, gtsv, gtui, gtvi,        & 
    88                                &                                   tsb, tsa, jpts        )  ! iso-level bilaplacian 
    89       CASE ( 3 )   ;   CALL tra_ldf_bilapg  ( kt, nit000, 'TRA',             tsb, tsa, jpts        )  ! s-coord. geopot. bilap. 
     72      ! 
     73      SELECT CASE ( nldf )                     !* compute lateral mixing trend and add it to the general trend 
     74      ! 
     75      CASE ( n_lap   )                                   ! laplacian: iso-level operator 
     76         CALL tra_ldf_lap  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb,      tsa, jpts,  1   ) 
    9077         ! 
    91       CASE ( -1 )                                ! esopa: test all possibility with control print 
    92          CALL tra_ldf_lap   ( kt, nit000, 'TRA', gtsu, gtsv, gtui, gtvi,        & 
    93          &                                       tsb, tsa, jpts        )  
    94          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf0 - Ta: ', mask1=tmask,               & 
    95          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    96          IF( ln_traldf_grif ) THEN 
    97             CALL tra_ldf_iso_grif( kt, nit000, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 ) 
    98          ELSE 
    99             CALL tra_ldf_iso     ( kt, nit000, 'TRA', gtsu, gtsv, gtui, gtvi,        & 
    100             &                                               tsb, tsa, jpts, ahtb0 )   
    101          ENDIF 
    102          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf1 - Ta: ', mask1=tmask,               & 
    103          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    104          CALL tra_ldf_bilap ( kt, nit000, 'TRA', gtsu, gtsv, gtui, gtvi,        & 
    105          &                                       tsb, tsa, jpts        )  
    106          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf2 - Ta: ', mask1=tmask,               & 
    107          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    108          CALL tra_ldf_bilapg( kt, nit000, 'TRA',             tsb, tsa, jpts        )  
    109          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf3 - Ta: ', mask1=tmask,               & 
    110          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
     78      CASE ( n_lap_i )                                   ! laplacian: standard iso-neutral operator (Madec) 
     79         CALL tra_ldf_iso  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts,  1   ) 
     80         ! 
     81      CASE ( n_lap_it )                                  ! laplacian: triad iso-neutral operator (griffies) 
     82         CALL tra_ldf_triad( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts,  1   ) 
     83         ! 
     84      CASE ( n_blp , n_blp_i , n_blp_it )                ! bilaplacian: iso-level & iso-neutral operators 
     85         CALL tra_ldf_blp  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb      , tsa, jpts, nldf ) 
    11186      END SELECT 
    11287 
    113 #if defined key_traldf_ano 
    114       tsa(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) - t0_ldf(:,:,:)      ! anomaly: substract the reference diffusivity 
    115       tsa(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) - s0_ldf(:,:,:) 
    116 #endif 
    117  
    118       IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics 
     88      IF( l_trdtra )   THEN                    !* save the horizontal diffusive trends for further diagnostics 
    11989         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    12090         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
    12191         CALL trd_tra( kt, 'TRA', jp_tem, jptra_ldf, ztrdt ) 
    12292         CALL trd_tra( kt, 'TRA', jp_sal, jptra_ldf, ztrds ) 
    123          CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )  
    124       ENDIF 
    125       !                                          ! print mean trends (used for debugging) 
     93         CALL wrk_dealloc( jpi,jpj,jpk,  ztrdt, ztrds )  
     94      ENDIF 
     95      !                                        !* print mean trends (used for debugging) 
    12696      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf  - Ta: ', mask1=tmask,               & 
    12797         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    12898      ! 
    129       IF( nn_timing == 1 )  CALL timing_stop('tra_ldf') 
     99      IF( nn_timing == 1 )   CALL timing_stop('tra_ldf') 
    130100      ! 
    131101   END SUBROUTINE tra_ldf 
    … …  
    139109      !! 
    140110      !! ** Method  :   set nldf from the namtra_ldf logicals 
    141       !!      nldf == -1   ESOPA test: ALL operators are used 
    142       !!      nldf ==  0   laplacian operator 
    143       !!      nldf ==  1   Rotated laplacian operator 
    144       !!      nldf ==  2   bilaplacian operator 
    145       !!      nldf ==  3   Rotated bilaplacian 
    146       !!---------------------------------------------------------------------- 
    147       INTEGER ::   ioptio, ierr         ! temporary integers  
    148       !!---------------------------------------------------------------------- 
    149  
    150       !  Define the lateral mixing oparator for tracers 
    151       ! =============================================== 
    152      
    153       IF(lwp) THEN                    ! Namelist print 
     111      !!---------------------------------------------------------------------- 
     112      INTEGER ::   ioptio, ierr   ! temporary integers  
     113      !!---------------------------------------------------------------------- 
     114      ! 
     115      IF(lwp) THEN                     ! Namelist print 
    154116         WRITE(numout,*) 
    155117         WRITE(numout,*) 'tra_ldf_init : lateral tracer diffusive operator' 
    … …  
    159121         WRITE(numout,*) 
    160122      ENDIF 
    161  
    162       !                               ! control the input 
     123      !                                ! control the input 
    163124      ioptio = 0 
    164       IF( ln_traldf_lap   )   ioptio = ioptio + 1 
    165       IF( ln_traldf_bilap )   ioptio = ioptio + 1 
    166       IF( ioptio >  1 )   CALL ctl_stop( '          use ONE or NONE of the 2 lap/bilap operator type on tracer' ) 
    167       IF( ioptio == 0 )   nldf = -2   ! No lateral diffusion 
     125      IF( ln_traldf_lap )   ioptio = ioptio + 1 
     126      IF( ln_traldf_blp )   ioptio = ioptio + 1 
     127      IF( ioptio >  1   )   CALL ctl_stop( 'tra_ldf_init: use ONE or NONE of the 2 lap/bilap operator type on tracer' ) 
     128      IF( ioptio == 0   )   nldf = n_no_ldf     ! No lateral diffusion 
    168129      ioptio = 0 
    169       IF( ln_traldf_level )   ioptio = ioptio + 1 
    170       IF( ln_traldf_hor   )   ioptio = ioptio + 1 
    171       IF( ln_traldf_iso   )   ioptio = ioptio + 1 
    172       IF( ioptio >  1 )   CALL ctl_stop( '          use only ONE direction (level/hor/iso)' ) 
    173  
    174       ! defined the type of lateral diffusion from ln_traldf_... logicals 
    175       ! CAUTION : nldf = 1 is used in trazdf_imp, change it carefully 
     130      IF( ln_traldf_lev )   ioptio = ioptio + 1 
     131      IF( ln_traldf_hor )   ioptio = ioptio + 1 
     132      IF( ln_traldf_iso )   ioptio = ioptio + 1 
     133      IF( ioptio >  1 )   CALL ctl_stop( 'tra_ldf_init: use only ONE direction (level/hor/iso)' ) 
     134      ! 
     135      !                                ! defined the type of lateral diffusion from ln_traldf_... logicals 
    176136      ierr = 0 
    177       IF( ln_traldf_lap ) THEN       ! laplacian operator 
    178          IF ( ln_zco ) THEN                ! z-coordinate 
    179             IF ( ln_traldf_level )   nldf = 0      ! iso-level  (no rotation) 
    180             IF ( ln_traldf_hor   )   nldf = 0      ! horizontal (no rotation) 
    181             IF ( ln_traldf_iso   )   nldf = 1      ! isoneutral (   rotation) 
    182          ENDIF 
    183          IF ( ln_zps ) THEN             ! zps-coordinate 
    184             IF ( ln_traldf_level )   ierr = 1      ! iso-level not allowed 
    185             IF ( ln_traldf_hor   )   nldf = 0      ! horizontal (no rotation) 
    186             IF ( ln_traldf_iso   )   nldf = 1      ! isoneutral (   rotation) 
    187          ENDIF 
    188          IF ( ln_sco ) THEN             ! s-coordinate 
    189             IF ( ln_traldf_level )   nldf = 0      ! iso-level  (no rotation) 
    190             IF ( ln_traldf_hor   )   nldf = 1      ! horizontal (   rotation) 
    191             IF ( ln_traldf_iso   )   nldf = 1      ! isoneutral (   rotation) 
    192          ENDIF 
    193       ENDIF 
    194  
    195       IF( ln_traldf_bilap ) THEN      ! bilaplacian operator 
    196          IF ( ln_zco ) THEN                ! z-coordinate 
    197             IF ( ln_traldf_level )   nldf = 2      ! iso-level  (no rotation) 
    198             IF ( ln_traldf_hor   )   nldf = 2      ! horizontal (no rotation) 
    199             IF ( ln_traldf_iso   )   ierr = 2      ! isoneutral (   rotation) 
    200          ENDIF 
    201          IF ( ln_zps ) THEN             ! zps-coordinate 
    202             IF ( ln_traldf_level )   ierr = 1      ! iso-level not allowed  
    203             IF ( ln_traldf_hor   )   nldf = 2      ! horizontal (no rotation) 
    204             IF ( ln_traldf_iso   )   ierr = 2      ! isoneutral (   rotation) 
    205          ENDIF 
    206          IF ( ln_sco ) THEN             ! s-coordinate 
    207             IF ( ln_traldf_level )   nldf = 2      ! iso-level  (no rotation) 
    208             IF ( ln_traldf_hor   )   nldf = 3      ! horizontal (   rotation) 
    209             IF ( ln_traldf_iso   )   ierr = 2      ! isoneutral (   rotation) 
    210          ENDIF 
    211       ENDIF 
    212  
    213       IF( nldf == 3 )   CALL ctl_warn( 'geopotential bilaplacian tracer diffusion in s-coords not thoroughly tested' ) 
     137      IF( ln_traldf_lap ) THEN         ! laplacian operator 
     138         IF ( ln_zco ) THEN               ! z-coordinate 
     139            IF ( ln_traldf_lev   )   nldf = n_lap     ! iso-level = horizontal (no rotation) 
     140            IF ( ln_traldf_hor   )   nldf = n_lap     ! iso-level = horizontal (no rotation) 
     141            IF ( ln_traldf_iso   )   nldf = n_lap_i   ! iso-neutral: standard  (   rotation) 
     142            IF ( ln_traldf_triad )   nldf = n_lap_it  ! iso-neutral: triad     (   rotation) 
     143         ENDIF 
     144         IF ( ln_zps ) THEN               ! z-coordinate with partial step 
     145            IF ( ln_traldf_lev   )   ierr = 1         ! iso-level not allowed  
     146            IF ( ln_traldf_hor   )   nldf = n_lap     ! horizontal (no rotation) 
     147            IF ( ln_traldf_iso   )   nldf = n_lap_i   ! iso-neutral: standard (rotation) 
     148            IF ( ln_traldf_triad )   nldf = n_lap_it  ! iso-neutral: triad    (rotation) 
     149         ENDIF 
     150         IF ( ln_sco ) THEN               ! s-coordinate 
     151            IF ( ln_traldf_lev   )   nldf = n_lap     ! iso-level  (no rotation) 
     152            IF ( ln_traldf_hor   )   nldf = n_lap_it  ! horizontal (   rotation)       !!gm   a checker.... 
     153            IF ( ln_traldf_iso   )   nldf = n_lap_i   ! iso-neutral: standard (rotation) 
     154            IF ( ln_traldf_triad )   nldf = n_lap_it  ! iso-neutral: triad    (rotation) 
     155         ENDIF 
     156      ENDIF 
     157      ! 
     158      IF( ln_traldf_blp ) THEN         ! bilaplacian operator 
     159         IF ( ln_zco ) THEN               ! z-coordinate 
     160            IF ( ln_traldf_lev   )   nldf = n_blp     ! iso-level = horizontal (no rotation) 
     161            IF ( ln_traldf_hor   )   nldf = n_blp     ! iso-level = horizontal (no rotation) 
     162            IF ( ln_traldf_iso   )   nldf = n_blp_i   ! iso-neutral: standard (rotation) 
     163            IF ( ln_traldf_triad )   nldf = n_blp_it  ! iso-neutral: triad    (rotation) 
     164         ENDIF 
     165         IF ( ln_zps ) THEN               ! z-coordinate with partial step 
     166            IF ( ln_traldf_lev   )   ierr = 1         ! iso-level not allowed  
     167            IF ( ln_traldf_hor   )   nldf = n_blp     ! horizontal (no rotation) 
     168            IF ( ln_traldf_iso   )   nldf = n_blp_i   ! iso-neutral: standard (rotation) 
     169            IF ( ln_traldf_triad )   nldf = n_blp_it  ! iso-neutral: triad    (rotation) 
     170         ENDIF 
     171         IF ( ln_sco ) THEN               ! s-coordinate 
     172            IF ( ln_traldf_lev   )   nldf = n_blp     ! iso-level  (no rotation) 
     173            IF ( ln_traldf_hor   )   nldf = n_blp_it  ! horizontal (   rotation)       !!gm   a checker.... 
     174            IF ( ln_traldf_iso   )   nldf = n_blp_i   ! iso-neutral: standard (rotation) 
     175            IF ( ln_traldf_triad )   nldf = n_blp_it  ! iso-neutral: triad    (rotation) 
     176         ENDIF 
     177      ENDIF 
     178      ! 
    214179      IF( ierr == 1 )   CALL ctl_stop( ' iso-level in z-coordinate - partial step, not allowed' ) 
    215       IF( ierr == 2 )   CALL ctl_stop( ' isoneutral bilaplacian operator does not exist' ) 
    216       IF( lk_traldf_eiv .AND. .NOT.ln_traldf_iso )   & 
    217            CALL ctl_stop( '          eddy induced velocity on tracers',   & 
    218            &              ' the eddy induced velocity on tracers requires isopycnal laplacian diffusion' ) 
    219       IF( nldf == 1 .OR. nldf == 3 ) THEN      ! rotation 
    220          IF( .NOT.lk_ldfslp )   CALL ctl_stop( '          the rotation of the diffusive tensor require key_ldfslp' ) 
    221          l_traldf_rot = .TRUE.                 ! needed for trazdf_imp 
    222       ENDIF 
    223  
    224       IF( lk_esopa ) THEN 
    225          IF(lwp) WRITE(numout,*) '          esopa control: use all lateral physics options' 
    226          nldf = -1 
    227       ENDIF 
    228  
     180      IF( ln_ldfeiv .AND. .NOT.( ln_traldf_iso .OR. ln_traldf_triad ) )                                    & 
     181           &            CALL ctl_stop( '          eddy induced velocity on tracers requires isopycnal',    & 
     182           &                                                                    ' laplacian diffusion' ) 
     183      IF(  nldf == n_lap_i .OR. nldf == n_lap_it .OR. & 
     184         & nldf == n_blp_i .OR. nldf == n_blp_it  )   l_ldfslp = .TRUE.    ! slope of neutral surfaces required  
     185      ! 
    229186      IF(lwp) THEN 
    230187         WRITE(numout,*) 
    231          IF( nldf == -2 )   WRITE(numout,*) '          NO lateral diffusion' 
    232          IF( nldf == -1 )   WRITE(numout,*) '          ESOPA test All scheme used' 
    233          IF( nldf ==  0 )   WRITE(numout,*) '          laplacian operator' 
    234          IF( nldf ==  1 )   WRITE(numout,*) '          Rotated laplacian operator' 
    235          IF( nldf ==  2 )   WRITE(numout,*) '          bilaplacian operator' 
    236          IF( nldf ==  3 )   WRITE(numout,*) '          Rotated bilaplacian' 
    237       ENDIF 
    238  
    239       ! Reference T & S diffusivity (if necessary) 
    240       ! =========================== 
    241       CALL ldf_ano 
     188         IF( nldf == n_no_ldf )   WRITE(numout,*) '          NO lateral diffusion' 
     189         IF( nldf == n_lap    )   WRITE(numout,*) '          laplacian iso-level operator' 
     190         IF( nldf == n_lap_i  )   WRITE(numout,*) '          Rotated laplacian operator (standard)' 
     191         IF( nldf == n_lap_it )   WRITE(numout,*) '          Rotated laplacian operator (triad)' 
     192         IF( nldf == n_blp    )   WRITE(numout,*) '          bilaplacian iso-level operator' 
     193         IF( nldf == n_blp_i  )   WRITE(numout,*) '          Rotated bilaplacian operator (standard)' 
     194         IF( nldf == n_blp_it )   WRITE(numout,*) '          Rotated bilaplacian operator (triad)' 
     195      ENDIF 
    242196      ! 
    243197   END SUBROUTINE tra_ldf_init 
    244  
    245 #if defined key_traldf_ano 
    246    !!---------------------------------------------------------------------- 
    247    !!   'key_traldf_ano'               T & S lateral diffusion on anomalies 
    248    !!---------------------------------------------------------------------- 
    249  
    250    SUBROUTINE ldf_ano 
    251       !!---------------------------------------------------------------------- 
    252       !!                  ***  ROUTINE ldf_ano  *** 
    253       !! 
    254       !! ** Purpose :   initializations of  
    255       !!---------------------------------------------------------------------- 
    256       ! 
    257       USE zdf_oce         ! vertical mixing 
    258       USE trazdf          ! vertical mixing: double diffusion 
    259       USE zdfddm          ! vertical mixing: double diffusion 
    260       ! 
    261       INTEGER  ::   jk              ! Dummy loop indice 
    262       INTEGER  ::   ierr            ! local integer 
    263       LOGICAL  ::   llsave          ! local logical 
    264       REAL(wp) ::   zt0, zs0, z12   ! local scalar 
    265       REAL(wp), POINTER, DIMENSION(:,:,:) :: zt_ref, zs_ref, ztb, zsb, zavt      
    266       !!---------------------------------------------------------------------- 
    267       ! 
    268       IF( nn_timing == 1 )  CALL timing_start('ldf_ano') 
    269       ! 
    270       CALL wrk_alloc( jpi, jpj, jpk, zt_ref, zs_ref, ztb, zsb, zavt )  
    271       ! 
    272  
    273       IF(lwp) THEN 
    274          WRITE(numout,*) 
    275          WRITE(numout,*) 'tra:ldf_ano : lateral diffusion acting on anomalies' 
    276          WRITE(numout,*) '~~~~~~~~~~~' 
    277       ENDIF 
    278  
    279       !                              ! allocate trabbl arrays 
    280       ALLOCATE( t0_ldf(jpi,jpj,jpk) , s0_ldf(jpi,jpj,jpk) , STAT=ierr ) 
    281       IF( lk_mpp    )   CALL mpp_sum( ierr ) 
    282       IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_ano: unable to allocate arrays' ) 
    283  
    284       ! defined the T & S reference profiles 
    285       ! ------------------------------------ 
    286       zt0 =10.e0                               ! homogeneous ocean 
    287       zs0 =35.e0 
    288       zt_ref(:,:,:) = 10.0 * tmask(:,:,:) 
    289       zs_ref(:,:,:) = 35.0 * tmask(:,:,:) 
    290       IF(lwp) WRITE(numout,*) '              homogeneous ocean T = ', zt0, ' S = ',zs0 
    291  
    292       ! Initialisation of gtui/gtvi in case of no cavity 
    293       IF ( .NOT. ln_isfcav ) THEN 
    294          gtui(:,:,:) = 0.0_wp 
    295          gtvi(:,:,:) = 0.0_wp 
    296       END IF 
    297       !                                        ! T & S profile (to be coded +namelist parameter 
    298  
    299       ! prepare the ldf computation 
    300       ! --------------------------- 
    301       llsave = l_trdtra 
    302       l_trdtra = .false.      ! desactivate trend computation 
    303       t0_ldf(:,:,:) = 0.e0 
    304       s0_ldf(:,:,:) = 0.e0 
    305       ztb   (:,:,:) = tsb (:,:,:,jp_tem) 
    306       zsb   (:,:,:) = tsb (:,:,:,jp_sal) 
    307       ua    (:,:,:) = tsa (:,:,:,jp_tem) 
    308       va    (:,:,:) = tsa (:,:,:,jp_sal) 
    309       zavt  (:,:,:) = avt(:,:,:) 
    310       IF( lk_zdfddm ) THEN CALL ctl_stop( ' key_traldf_ano with key_zdfddm not implemented' ) 
    311       ! set tb, sb to reference values and avr to zero 
    312       tsb (:,:,:,jp_tem) = zt_ref(:,:,:) 
    313       tsb (:,:,:,jp_sal) = zs_ref(:,:,:) 
    314       tsa (:,:,:,jp_tem) = 0.e0 
    315       tsa (:,:,:,jp_sal) = 0.e0 
    316       avt(:,:,:)         = 0.e0 
    317  
    318       ! Compute the ldf trends 
    319       ! ---------------------- 
    320       CALL tra_ldf( nit000 + 1 )      ! horizontal components (+1: no more init) 
    321       CALL tra_zdf( nit000     )      ! vertical component (if necessary nit000 to performed the init) 
    322  
    323       ! finalise the computation and recover all arrays 
    324       ! ----------------------------------------------- 
    325       l_trdtra = llsave 
    326       z12 = 2.e0 
    327       IF( neuler == 1)   z12 = 1.e0 
    328       IF( ln_zdfexp ) THEN      ! ta,sa are the trends 
    329          t0_ldf(:,:,:) = tsa(:,:,:,jp_tem) 
    330          s0_ldf(:,:,:) = tsa(:,:,:,jp_sal) 
    331       ELSE 
    332          DO jk = 1, jpkm1 
    333             t0_ldf(:,:,jk) = ( tsa(:,:,jk,jp_tem) - tsb(:,:,jk,jp_tem) ) / ( z12 *rdttra(jk) ) 
    334             s0_ldf(:,:,jk) = ( tsa(:,:,jk,jp_sal) - tsb(:,:,jk,jp_sal) ) / ( z12 *rdttra(jk) ) 
    335          END DO 
    336       ENDIF 
    337       tsb(:,:,:,jp_tem) = ztb (:,:,:) 
    338       tsb(:,:,:,jp_sal) = zsb (:,:,:) 
    339       tsa(:,:,:,jp_tem) = ua  (:,:,:) 
    340       tsa(:,:,:,jp_sal) = va  (:,:,:) 
    341       avt(:,:,:)        = zavt(:,:,:) 
    342       ! 
    343       CALL wrk_dealloc( jpi, jpj, jpk, zt_ref, zs_ref, ztb, zsb, zavt )  
    344       ! 
    345       IF( nn_timing == 1 )  CALL timing_stop('ldf_ano') 
    346       ! 
    347    END SUBROUTINE ldf_ano 
    348  
    349 #else 
    350    !!---------------------------------------------------------------------- 
    351    !!   default option :   Dummy code   NO T & S background profiles 
    352    !!---------------------------------------------------------------------- 
    353    SUBROUTINE ldf_ano 
    354       IF(lwp) THEN 
    355          WRITE(numout,*) 
    356          WRITE(numout,*) 'tra:ldf_ano : lateral diffusion acting on the full fields' 
    357          WRITE(numout,*) '~~~~~~~~~~~' 
    358       ENDIF 
    359    END SUBROUTINE ldf_ano 
    360 #endif 
    361198 
    362199   !!====================================================================== 

  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90

    r5737 r5758  
    44   !! Ocean  tracers:  horizontal component of the lateral tracer mixing trend 
    55   !!====================================================================== 
    6    !! History :  OPA  !  1994-08  (G. Madec, M. Imbard) 
    7    !!            8.0  !  1997-05  (G. Madec)  split into traldf and trazdf 
    8    !!            NEMO !  2002-08  (G. Madec)  Free form, F90 
    9    !!            1.0  !  2005-11  (G. Madec)  merge traldf and trazdf :-) 
    10    !!            3.3  !  2010-09  (C. Ethe, G. Madec) Merge TRA-TRC 
     6   !! History :  OPA  ! 1994-08  (G. Madec, M. Imbard) 
     7   !!            8.0  ! 1997-05  (G. Madec)  split into traldf and trazdf 
     8   !!            NEMO ! 2002-08  (G. Madec)  Free form, F90 
     9   !!            1.0  ! 2005-11  (G. Madec)  merge traldf and trazdf :-) 
     10   !!            3.3  ! 2010-09  (C. Ethe, G. Madec) Merge TRA-TRC 
     11   !!            3.7  ! 2014-01  (G. Madec, S. Masson)  restructuration/simplification of aht/aeiv specification 
     12   !!             -   ! 2014-02  (F. Lemarie, G. Madec)  triad operator (Griffies) + Method of Stabilizing Correction 
    1113   !!---------------------------------------------------------------------- 
    12 #if   defined key_ldfslp   ||   defined key_esopa 
     14 
    1315   !!---------------------------------------------------------------------- 
    14    !!   'key_ldfslp'               slope of the lateral diffusive direction 
    15    !!---------------------------------------------------------------------- 
    16    !!   tra_ldf_iso  : update the tracer trend with the horizontal  
    17    !!                  component of a iso-neutral laplacian operator 
    18    !!                  and with the vertical part of 
    19    !!                  the isopycnal or geopotential s-coord. operator  
     16   !!   tra_ldf_iso  : update the tracer trend with the horizontal component of a iso-neutral laplacian operator 
     17   !!                  and with the vertical part of the isopycnal or geopotential s-coord. operator  
    2018   !!---------------------------------------------------------------------- 
    2119   USE oce             ! ocean dynamics and active tracers 
    … …  
    2321   USE trc_oce         ! share passive tracers/Ocean variables 
    2422   USE zdf_oce         ! ocean vertical physics 
    25    USE ldftra_oce      ! ocean active tracers: lateral physics 
     23   USE ldftra          ! lateral diffusion: tracer eddy coefficients 
    2624   USE ldfslp          ! iso-neutral slopes 
    2725   USE diaptr          ! poleward transport diagnostics 
     26   ! 
    2827   USE in_out_manager  ! I/O manager 
    2928   USE iom             ! I/O library 
    … …  
    4039   !! * Substitutions 
    4140#  include "domzgr_substitute.h90" 
    42 #  include "ldftra_substitute.h90" 
    4341#  include "vectopt_loop_substitute.h90" 
    4442   !!---------------------------------------------------------------------- 
    45    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     43   !! NEMO/OPA 3.7 , NEMO Consortium (2015) 
    4644   !! $Id$ 
    4745   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    4947CONTAINS 
    5048 
    51    SUBROUTINE tra_ldf_iso( kt, kit000, cdtype, pgu, pgv,              & 
    52       &                                pgui, pgvi,                    & 
    53       &                                ptb, pta, kjpt, pahtb0 ) 
     49  SUBROUTINE tra_ldf_iso( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   & 
     50      &                                                   pgui, pgvi,   & 
     51      &                                       ptb , ptbb, pta , kjpt, kpass ) 
    5452      !!---------------------------------------------------------------------- 
    5553      !!                  ***  ROUTINE tra_ldf_iso  *** 
    … …  
    6664      !! 
    6765      !!      1st part :  masked horizontal derivative of T  ( di[ t ] ) 
    68       !!      ========    with partial cell update if ln_zps=T. 
     66      !!      ========    with partial cell update if ln_zps=T 
     67      !!                  with top     cell update if ln_isfcav 
    6968      !! 
    7069      !!      2nd part :  horizontal fluxes of the lateral mixing operator 
    7170      !!      ========     
    72       !!         zftu = (aht+ahtb0) e2u*e3u/e1u di[ tb ] 
    73       !!               - aht      e2u*uslp    dk[ mi(mk(tb)) ] 
    74       !!         zftv = (aht+ahtb0) e1v*e3v/e2v dj[ tb ] 
    75       !!               - aht      e2u*vslp    dk[ mj(mk(tb)) ] 
     71      !!         zftu =  pahu e2u*e3u/e1u di[ tb ] 
     72      !!               - pahu e2u*uslp    dk[ mi(mk(tb)) ] 
     73      !!         zftv =  pahv e1v*e3v/e2v dj[ tb ] 
     74      !!               - pahv e2u*vslp    dk[ mj(mk(tb)) ] 
    7675      !!      take the horizontal divergence of the fluxes: 
    77       !!         difft = 1/(e1t*e2t*e3t) {  di-1[ zftu ] +  dj-1[ zftv ]  } 
     76      !!         difft = 1/(e1e2t*e3t) {  di-1[ zftu ] +  dj-1[ zftv ]  } 
    7877      !!      Add this trend to the general trend (ta,sa): 
    7978      !!         ta = ta + difft 
    … …  
    8281      !!      ========  (excluding the vertical flux proportional to dk[t] ) 
    8382      !!      vertical fluxes associated with the rotated lateral mixing: 
    84       !!         zftw =-aht { e2t*wslpi di[ mi(mk(tb)) ] 
    85       !!                     + e1t*wslpj dj[ mj(mk(tb)) ]  } 
     83      !!         zftw = - {  mi(mk(pahu)) * e2t*wslpi di[ mi(mk(tb)) ] 
     84      !!                   + mj(mk(pahv)) * e1t*wslpj dj[ mj(mk(tb)) ]  } 
    8685      !!      take the horizontal divergence of the fluxes: 
    87       !!         difft = 1/(e1t*e2t*e3t) dk[ zftw ] 
     86      !!         difft = 1/(e1e2t*e3t) dk[ zftw ] 
    8887      !!      Add this trend to the general trend (ta,sa): 
    8988      !!         pta = pta + difft 
    … …  
    9190      !! ** Action :   Update pta arrays with the before rotated diffusion 
    9291      !!---------------------------------------------------------------------- 
    93       USE oce     , ONLY:   zftu => ua       , zftv  => va         ! (ua,va) used as workspace 
    94       ! 
    9592      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index 
    96       INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index 
     93      INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index 
    9794      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
    9895      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers 
    99       REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu , pgv    ! tracer gradient at pstep levels 
    100       REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgui, pgvi   ! tracer gradient at pstep levels 
    101       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields 
    102       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend  
    103       REAL(wp)                             , INTENT(in   ) ::   pahtb0     ! background diffusion coef 
     96      INTEGER                              , INTENT(in   ) ::   kpass      ! =1/2 first or second passage 
     97      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s] 
     98      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels 
     99      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels 
     100      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! tracer (kpass=1) or laplacian of tracer (kpass=2) 
     101      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptbb       ! tracer (only used in kpass=2) 
     102      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend 
    104103      ! 
    105104      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
    106       INTEGER  ::  ikt 
    107       REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3   ! local scalars 
    108       REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4   !   -      - 
    109       REAL(wp) ::  zcoef0, zbtr, ztra            !   -      - 
    110       REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d 
    111       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zdkt, zdk1t, zdit, zdjt, ztfw  
     105      INTEGER  ::  ierr             ! local integer 
     106      REAL(wp) ::  zmsku, zahu_w, zabe1, zcof1, zcoef3   ! local scalars 
     107      REAL(wp) ::  zmskv, zahv_w, zabe2, zcof2, zcoef4   !   -      - 
     108      REAL(wp) ::  zcoef0, ze3w_2, zsign, z2dt, z1_2dt   !   -      - 
     109#if defined key_diaar5 
     110      REAL(wp) ::   zztmp   ! local scalar 
     111#endif 
     112      REAL(wp), POINTER, DIMENSION(:,:)   ::   zdkt, zdk1t, z2d 
     113      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdit, zdjt, zftu, zftv, ztfw  
    112114      !!---------------------------------------------------------------------- 
    113115      ! 
    114116      IF( nn_timing == 1 )  CALL timing_start('tra_ldf_iso') 
    115117      ! 
    116       CALL wrk_alloc( jpi, jpj,      z2d )  
    117       CALL wrk_alloc( jpi, jpj, jpk, zdit, zdjt, ztfw, zdkt, zdk1t )  
    118       ! 
    119  
     118      CALL wrk_alloc( jpi,jpj,       zdkt, zdk1t, z2d )  
     119      CALL wrk_alloc( jpi,jpj,jpk,   zdit, zdjt , zftu, zftv, ztfw  )  
     120      ! 
    120121      IF( kt == kit000 )  THEN 
    121122         IF(lwp) WRITE(numout,*) 
    122123         IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator on ', cdtype 
    123124         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
     125         ! 
     126         akz     (:,:,:) = 0._wp       
     127         ah_wslp2(:,:,:) = 0._wp 
     128      ENDIF 
     129      ! 
     130      !                                               ! set time step size (Euler/Leapfrog) 
     131      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   z2dt =     rdttra(1)      ! at nit000   (Euler) 
     132      ELSE                                        ;   z2dt = 2.* rdttra(1)      !             (Leapfrog) 
     133      ENDIF 
     134      z1_2dt = 1._wp / z2dt 
     135      ! 
     136      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0) 
     137      ELSE                    ;   zsign = -1._wp 
     138      ENDIF 
     139          
     140          
     141      !!---------------------------------------------------------------------- 
     142      !!   0 - calculate  ah_wslp2 and akz 
     143      !!---------------------------------------------------------------------- 
     144      ! 
     145      IF( kpass == 1 ) THEN                  !==  first pass only  ==! 
     146         ! 
     147         DO jk = 2, jpkm1 
     148            DO jj = 2, jpjm1 
     149               DO ji = fs_2, fs_jpim1   ! vector opt. 
     150                  ! 
     151                  zmsku = tmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          & 
     152                     &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  ) 
     153                  zmskv = tmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          & 
     154                     &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  ) 
     155                     ! 
     156                  zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    & 
     157                     &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku 
     158                  zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    & 
     159                     &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv 
     160                     ! 
     161                  ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   & 
     162                     &               + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) 
     163               END DO 
     164            END DO 
     165         END DO 
     166         ! 
     167         IF( ln_traldf_msc ) THEN                ! stabilizing vertical diffusivity coefficient 
     168            DO jk = 2, jpkm1 
     169               DO jj = 2, jpjm1 
     170                  DO ji = fs_2, fs_jpim1 
     171                     akz(ji,jj,jk) = 0.25_wp * (                                                                     & 
     172                        &              ( pahu(ji  ,jj,jk) + pahu(ji  ,jj,jk-1) ) / ( e1u(ji  ,jj) * e1u(ji  ,jj) )   & 
     173                        &            + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) )   & 
     174                        &            + ( pahv(ji,jj  ,jk) + pahv(ji,jj  ,jk-1) ) / ( e2v(ji,jj  ) * e2v(ji,jj  ) )   & 
     175                        &            + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) )   ) 
     176                  END DO 
     177               END DO 
     178            END DO 
     179            ! 
     180            IF( ln_traldf_blp ) THEN                ! bilaplacian operator 
     181               DO jk = 2, jpkm1 
     182                  DO jj = 1, jpjm1 
     183                     DO ji = 1, fs_jpim1 
     184                        akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   & 
     185                           &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( fse3w(ji,jj,jk) * fse3w(ji,jj,jk) )  ) 
     186                     END DO 
     187                  END DO 
     188               END DO 
     189            ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator 
     190               DO jk = 2, jpkm1 
     191                  DO jj = 1, jpjm1 
     192                     DO ji = 1, fs_jpim1 
     193                        ze3w_2 = fse3w(ji,jj,jk) * fse3w(ji,jj,jk) 
     194                        zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
     195                        akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt 
     196                     END DO 
     197                  END DO 
     198               END DO 
     199           ENDIF 
     200           ! 
     201         ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit 
     202            akz(:,:,:) = ah_wslp2(:,:,:)       
     203         ENDIF 
    124204      ENDIF 
    125205      ! 
    … …  
    131211         !!   I - masked horizontal derivative  
    132212         !!---------------------------------------------------------------------- 
    133          !!bug ajout.... why?   ( 1,jpj,:) and (jpi,1,:) should be sufficient.... 
    134          zdit (1,:,:) = 0.e0     ;     zdit (jpi,:,:) = 0.e0 
    135          zdjt (1,:,:) = 0.e0     ;     zdjt (jpi,:,:) = 0.e0 
     213!!gm : bug.... why (x,:,:)?   (1,jpj,:) and (jpi,1,:) should be sufficient.... 
     214         zdit (1,:,:) = 0._wp     ;     zdit (jpi,:,:) = 0._wp 
     215         zdjt (1,:,:) = 0._wp     ;     zdjt (jpi,:,:) = 0._wp 
    136216         !!end 
    137217 
    … …  
    145225            END DO 
    146226         END DO 
    147  
    148          ! partial cell correction 
    149          IF( ln_zps ) THEN      ! partial steps correction at the last ocean level  
    150             DO jj = 1, jpjm1 
     227         IF( ln_zps ) THEN      ! botton and surface ocean correction of the horizontal gradient 
     228            DO jj = 1, jpjm1              ! bottom correction (partial bottom cell) 
    151229               DO ji = 1, fs_jpim1   ! vector opt. 
    152 ! IF useless if zpshde defines pgu everywhere 
     230!!gm  the following anonymous remark is to considered: ! IF useless if zpshde defines pgu everywhere 
    153231                  zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)           
    154232                  zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 
    155233               END DO 
    156234            END DO 
     235            IF( ln_isfcav ) THEN      ! first wet level beneath a cavity 
     236               DO jj = 1, jpjm1 
     237                  DO ji = 1, fs_jpim1   ! vector opt. 
     238                     IF( miku(ji,jj) > 1 )   zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)           
     239                     IF( mikv(ji,jj) > 1 )   zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)      
     240                  END DO 
     241               END DO 
     242            ENDIF 
    157243         ENDIF 
    158          IF( ln_zps .AND. ln_isfcav ) THEN      ! partial steps correction at the first wet level beneath a cavity 
    159             DO jj = 1, jpjm1 
     244 
     245         !!---------------------------------------------------------------------- 
     246         !!   II - horizontal trend  (full) 
     247         !!---------------------------------------------------------------------- 
     248         ! 
     249         DO jk = 1, jpkm1                                 ! Horizontal slab 
     250            ! 
     251            !                             !== Vertical tracer gradient 
     252            zdk1t(:,:) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * wmask(:,:,jk+1)     ! level jk+1 
     253            ! 
     254            IF( jk == 1 ) THEN   ;   zdkt(:,:) = zdk1t(:,:)                          ! surface: zdkt(jk=1)=zdkt(jk=2) 
     255            ELSE                 ;   zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk) 
     256            ENDIF 
     257!!gm I don't understand why we should need this.... since wmask is used instead of tmask 
     258!         IF ( ln_isfcav ) THEN 
     259!            DO jj = 1, jpj 
     260!               DO ji = 1, jpi   ! vector opt. 
     261!                  ikt = mikt(ji,jj) ! surface level 
     262!                  zdk1t(ji,jj,ikt) = ( ptb(ji,jj,ikt,jn  ) - ptb(ji,jj,ikt+1,jn) ) * wmask(ji,jj,ikt+1) 
     263!                  zdkt (ji,jj,ikt) = zdk1t(ji,jj,ikt) 
     264!               END DO 
     265!            END DO 
     266!         END IF 
     267!!gm  
     268 
     269            DO jj = 1 , jpjm1            !==  Horizontal fluxes 
    160270               DO ji = 1, fs_jpim1   ! vector opt. 
    161                   IF (miku(ji,jj) > 1) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)           
    162                   IF (mikv(ji,jj) > 1) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)      
    163                END DO 
    164             END DO 
    165          END IF 
    166  
    167          !!---------------------------------------------------------------------- 
    168          !!   II - horizontal trend  (full) 
    169          !!---------------------------------------------------------------------- 
    170 !!!!!!!!!!CDIR PARALLEL DO PRIVATE( zdk1t )  
    171             ! 1. Vertical tracer gradient at level jk and jk+1 
    172             ! ------------------------------------------------ 
    173          !  
    174          ! interior value  
    175          DO jk = 2, jpkm1                
    176             DO jj = 1, jpj 
    177                DO ji = 1, jpi   ! vector opt. 
    178                   zdk1t(ji,jj,jk) = ( ptb(ji,jj,jk,jn  ) - ptb(ji,jj,jk+1,jn) ) * wmask(ji,jj,jk+1) 
    179                   ! 
    180                   zdkt(ji,jj,jk)  = ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn  ) ) * wmask(ji,jj,jk) 
    181                END DO 
    182             END DO 
    183          END DO 
    184          ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 
    185          zdk1t(:,:,1) = ( ptb(:,:,1,jn  ) - ptb(:,:,2,jn) ) * wmask(:,:,2) 
    186          zdkt (:,:,1) = zdk1t(:,:,1) 
    187          IF ( ln_isfcav ) THEN 
    188             DO jj = 1, jpj 
    189                DO ji = 1, jpi   ! vector opt. 
    190                   ikt = mikt(ji,jj) ! surface level 
    191                   zdk1t(ji,jj,ikt) = ( ptb(ji,jj,ikt,jn  ) - ptb(ji,jj,ikt+1,jn) ) * wmask(ji,jj,ikt+1) 
    192                   zdkt (ji,jj,ikt) = zdk1t(ji,jj,ikt) 
    193                END DO 
    194             END DO 
    195          END IF 
    196  
    197          ! 2. Horizontal fluxes 
    198          ! --------------------    
    199          DO jk = 1, jpkm1 
    200             DO jj = 1 , jpjm1 
    201                DO ji = 1, fs_jpim1   ! vector opt. 
    202                   zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * e2_e1u(ji,jj) * fse3u_n(ji,jj,jk) 
    203                   zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * e1_e2v(ji,jj) * fse3v_n(ji,jj,jk) 
     271                  zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * fse3u_n(ji,jj,jk) 
     272                  zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * fse3v_n(ji,jj,jk) 
    204273                  ! 
    205274                  zmsku = 1. / MAX(  tmask(ji+1,jj,jk  ) + tmask(ji,jj,jk+1)   & 
    … …  
    209278                     &             + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk  ), 1. ) 
    210279                  ! 
    211                   zcof1 = - fsahtu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku 
    212                   zcof2 = - fsahtv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv 
     280                  zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku 
     281                  zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv 
    213282                  ! 
    214283                  zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)   & 
    215                      &              + zcof1 * (  zdkt (ji+1,jj,jk) + zdk1t(ji,jj,jk)      & 
    216                      &                         + zdk1t(ji+1,jj,jk) + zdkt (ji,jj,jk)  )  ) * umask(ji,jj,jk) 
     284                     &               + zcof1 * (  zdkt (ji+1,jj) + zdk1t(ji,jj)      & 
     285                     &                          + zdk1t(ji+1,jj) + zdkt (ji,jj)  )  ) * umask(ji,jj,jk) 
    217286                  zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   & 
    218                      &              + zcof2 * (  zdkt (ji,jj+1,jk) + zdk1t(ji,jj,jk)      & 
    219                      &                         + zdk1t(ji,jj+1,jk) + zdkt (ji,jj,jk)  )  ) * vmask(ji,jj,jk)                   
    220                END DO 
    221             END DO 
    222  
    223             ! II.4 Second derivative (divergence) and add to the general trend 
    224             ! ---------------------------------------------------------------- 
    225             DO jj = 2 , jpjm1 
     287                     &               + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj)      & 
     288                     &                          + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk)                   
     289               END DO 
     290            END DO 
     291            ! 
     292            DO jj = 2 , jpjm1          !== horizontal divergence and add to pta 
    226293               DO ji = fs_2, fs_jpim1   ! vector opt. 
    227                   zbtr = 1.0 / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 
    228                   ztra = zbtr * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  ) 
    229                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 
    230                END DO 
    231             END DO 
    232             !                                          ! =============== 
     294                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      & 
     295                     &                                           + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   & 
     296                     &                                        / (  e1e2t(ji,jj) * fse3t(ji,jj,jk)  ) 
     297               END DO 
     298            END DO 
    233299         END DO                                        !   End of slab   
    234          !                                             ! =============== 
    235          ! 
    236          ! "Poleward" diffusive heat or salt transports (T-S case only) 
    237          IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
    238             ! note sign is reversed to give down-gradient diffusive transports (#1043) 
    239             IF( jn == jp_tem)   htr_ldf(:) = ptr_sj( -zftv(:,:,:) ) 
    240             IF( jn == jp_sal)   str_ldf(:) = ptr_sj( -zftv(:,:,:) ) 
    241          ENDIF 
    242   
    243          IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN 
    244            ! 
    245            IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN 
    246                z2d(:,:) = 0._wp  
    247                DO jk = 1, jpkm1 
    248                   DO jj = 2, jpjm1 
    249                      DO ji = fs_2, fs_jpim1   ! vector opt. 
    250                         z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk)  
    251                      END DO 
    252                   END DO 
    253                END DO 
    254                z2d(:,:) = - rau0_rcp * z2d(:,:)     ! note sign is reversed to give down-gradient diffusive transports (#1043) 
    255                CALL lbc_lnk( z2d, 'U', -1. ) 
    256                CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction 
    257                ! 
    258                z2d(:,:) = 0._wp  
    259                DO jk = 1, jpkm1 
    260                   DO jj = 2, jpjm1 
    261                      DO ji = fs_2, fs_jpim1   ! vector opt. 
    262                         z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk)  
    263                      END DO 
    264                   END DO 
    265                END DO 
    266                z2d(:,:) = - rau0_rcp * z2d(:,:)     ! note sign is reversed to give down-gradient diffusive transports (#1043) 
    267                CALL lbc_lnk( z2d, 'V', -1. ) 
    268                CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in i-direction 
    269             END IF 
    270             ! 
    271          ENDIF 
    272  
    273          !!---------------------------------------------------------------------- 
    274          !!   III - vertical trend of T & S (extra diagonal terms only) 
    275          !!---------------------------------------------------------------------- 
    276           
    277          ! Local constant initialization 
    278          ! ----------------------------- 
    279          ztfw(1,:,:) = 0.e0     ;     ztfw(jpi,:,:) = 0.e0 
     300 
     301 
     302         !!---------------------------------------------------------------------- 
     303         !!   III - vertical trend (full) 
     304         !!---------------------------------------------------------------------- 
     305          
     306         ztfw(1,:,:) = 0._wp     ;     ztfw(jpi,:,:) = 0._wp 
    280307          
    281308         ! Vertical fluxes 
    … …  
    283310          
    284311         ! Surface and bottom vertical fluxes set to zero 
    285          ztfw(:,:, 1 ) = 0.e0      ;      ztfw(:,:,jpk) = 0.e0 
     312         ztfw(:,:, 1 ) = 0._wp      ;      ztfw(:,:,jpk) = 0._wp 
    286313          
    287314         ! interior (2=<jk=<jpk-1) 
    … …  
    289316            DO jj = 2, jpjm1 
    290317               DO ji = fs_2, fs_jpim1   ! vector opt. 
    291                   zcoef0 = - fsahtw(ji,jj,jk) * wmask(ji,jj,jk) 
    292                   ! 
    293                   zmsku = 1./MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)      & 
    294                      &            + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk), 1.  ) 
    295                   zmskv = 1./MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)      & 
    296                      &            + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk), 1.  ) 
    297                   ! 
    298                   zcoef3 = zcoef0 * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk) 
    299                   zcoef4 = zcoef0 * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) 
     318                  ! 
     319                  zmsku = tmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          & 
     320                     &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  ) 
     321                  zmskv = tmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          & 
     322                     &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  ) 
     323                     ! 
     324                  zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    & 
     325                     &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku 
     326                  zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    & 
     327                     &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv 
     328                     ! 
     329                  zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk)   !wslpi & j are already w-masked 
     330                  zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) 
    300331                  ! 
    301332                  ztfw(ji,jj,jk) = zcoef3 * (   zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)      & 
    … …  
    306337            END DO 
    307338         END DO 
    308           
    309           
    310          ! I.5 Divergence of vertical fluxes added to the general tracer trend 
    311          ! ------------------------------------------------------------------- 
    312          DO jk = 1, jpkm1 
     339         ! 
     340         !                                !==  add the vertical 33 flux  ==! 
     341         IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz 
     342            DO jk = 2, jpkm1        
     343               DO jj = 1, jpjm1 
     344                  DO ji = fs_2, fs_jpim1 
     345                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk)   & 
     346                        &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             & 
     347                        &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 
     348                  END DO 
     349               END DO 
     350            END DO 
     351            ! 
     352         ELSE                                   ! bilaplacian  
     353            SELECT CASE( kpass ) 
     354            CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2 
     355               DO jk = 2, jpkm1  
     356                  DO jj = 1, jpjm1 
     357                     DO ji = fs_2, fs_jpim1 
     358                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk)    & 
     359                           &           + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   & 
     360                           &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * tmask(ji,jj,jk) / fse3w(ji,jj,jk) 
     361                     END DO 
     362                  END DO 
     363               END DO  
     364            CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb  and ptbb gradients, resp. 
     365               DO jk = 2, jpkm1  
     366                  DO jj = 1, jpjm1 
     367                     DO ji = fs_2, fs_jpim1 
     368                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk)                      & 
     369                           &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   & 
     370                           &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   ) 
     371                     END DO 
     372                  END DO 
     373               END DO 
     374            END SELECT 
     375         ENDIF 
     376         !          
     377         DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pta  ==! 
    313378            DO jj = 2, jpjm1 
    314379               DO ji = fs_2, fs_jpim1   ! vector opt. 
    315                   zbtr = 1.0 / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 
    316                   ztra = (  ztfw(ji,jj,jk) - ztfw(ji,jj,jk+1)  ) * zbtr 
    317                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 
     380                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   & 
     381                     &                                        / (  e1e2t(ji,jj) * fse3t_n(ji,jj,jk)  ) 
    318382               END DO 
    319383            END DO 
    320384         END DO 
    321385         ! 
    322       END DO 
    323       ! 
    324       CALL wrk_dealloc( jpi, jpj, z2d )  
    325       CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt, ztfw, zdkt, zdk1t )  
     386         IF( ( kpass == 1 .AND. ln_traldf_lap ) .OR.  &     !==  first pass only (  laplacian)  ==! 
     387             ( kpass == 2 .AND. ln_traldf_blp ) ) THEN      !==  2nd   pass      (bilaplacian)  ==! 
     388            ! 
     389            !                             ! "Poleward" diffusive heat or salt transports (T-S case only) 
     390            IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
     391               ! note sign is reversed to give down-gradient diffusive transports (#1043) 
     392               IF( jn == jp_tem)   htr_ldf(:) = ptr_sj( -zftv(:,:,:) ) 
     393               IF( jn == jp_sal)   str_ldf(:) = ptr_sj( -zftv(:,:,:) ) 
     394            ENDIF 
     395            ! 
     396            IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN 
     397              ! 
     398              IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN 
     399                  z2d(:,:) = zftu(ji,jj,1)  
     400                  DO jk = 2, jpkm1 
     401                     DO jj = 2, jpjm1 
     402                        DO ji = fs_2, fs_jpim1   ! vector opt. 
     403                           z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk)  
     404                        END DO 
     405                     END DO 
     406                  END DO 
     407!!gm CAUTION I think there is an error of sign when using BLP operator.... 
     408!!gm         a multiplication by zsign is required (to be checked twice !) 
     409                  z2d(:,:) = - rau0_rcp * z2d(:,:)     ! note sign is reversed to give down-gradient diffusive transports (#1043) 
     410                  CALL lbc_lnk( z2d, 'U', -1. ) 
     411                  CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction 
     412                  ! 
     413                  z2d(:,:) = zftv(ji,jj,1)  
     414                  DO jk = 2, jpkm1 
     415                     DO jj = 2, jpjm1 
     416                        DO ji = fs_2, fs_jpim1   ! vector opt. 
     417                           z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk)  
     418                        END DO 
     419                     END DO 
     420                  END DO 
     421                  z2d(:,:) = - rau0_rcp * z2d(:,:)     ! note sign is reversed to give down-gradient diffusive transports (#1043) 
     422                  CALL lbc_lnk( z2d, 'V', -1. ) 
     423                  CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in i-direction 
     424               END IF 
     425               ! 
     426            ENDIF 
     427            ! 
     428         ENDIF                                                    !== end pass selection  ==! 
     429         ! 
     430         !                                                        ! =============== 
     431      END DO                                                      ! end tracer loop 
     432      !                                                           ! =============== 
     433      ! 
     434      CALL wrk_dealloc( jpi, jpj,      zdkt, zdk1t, z2d )  
     435      CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt , zftu, zftv, ztfw  )  
    326436      ! 
    327437      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_iso') 
    328438      ! 
    329439   END SUBROUTINE tra_ldf_iso 
    330  
    331 #else 
    332    !!---------------------------------------------------------------------- 
    333    !!   default option :   Dummy code   NO rotation of the diffusive tensor 
    334    !!---------------------------------------------------------------------- 
    335 CONTAINS 
    336    SUBROUTINE tra_ldf_iso( kt, kit000,cdtype, pgu, pgv, pgui, pgvi, ptb, pta, kjpt, pahtb0 )      ! Empty routine 
    337       INTEGER:: kt, kit000 
    338       CHARACTER(len=3) ::   cdtype 
    339       REAL, DIMENSION(:,:,:) ::   pgu, pgv, pgui, pgvi    ! tracer gradient at pstep levels 
    340       REAL, DIMENSION(:,:,:,:) ::   ptb, pta 
    341       WRITE(*,*) 'tra_ldf_iso: You should not have seen this print! error?', kt, kit000, cdtype,   & 
    342          &                       pgu(1,1,1), pgv(1,1,1), ptb(1,1,1,1), pta(1,1,1,1), kjpt, pahtb0 
    343    END SUBROUTINE tra_ldf_iso 
    344 #endif 
    345440 
    346441   !!============================================================================== 

  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90

    r5737 r5758  
    22   !!============================================================================== 
    33   !!                       ***  MODULE  traldf_lap  *** 
    4    !! Ocean  tracers:  horizontal component of the lateral tracer mixing trend 
     4   !! Ocean tracers:  lateral diffusivity trend  (laplacian and bilaplacian) 
    55   !!============================================================================== 
    6    !! History :  OPA  !  87-06  (P. Andrich, D. L Hostis)  Original code 
    7    !!                 !  91-11  (G. Madec) 
    8    !!                 !  95-11  (G. Madec)  suppress volumetric scale factors 
    9    !!                 !  96-01  (G. Madec)  statement function for e3 
    10    !!            NEMO !  02-06  (G. Madec)  F90: Free form and module 
    11    !!            1.0  !  04-08  (C. Talandier) New trends organization 
    12    !!                 !  05-11  (G. Madec)  add zps case 
    13    !!            3.0  !  10-06  (C. Ethe, G. Madec) Merge TRA-TRC 
     6   !! History :  OPA  ! 1987-06  (P. Andrich, D. L Hostis)  Original code 
     7   !!                 ! 1991-11  (G. Madec) 
     8   !!                 ! 1995-11  (G. Madec)  suppress volumetric scale factors 
     9   !!                 ! 1996-01  (G. Madec)  statement function for e3 
     10   !!            NEMO ! 2002-06  (G. Madec)  F90: Free form and module 
     11   !!            1.0  ! 2004-08  (C. Talandier) New trends organization 
     12   !!                 ! 2005-11  (G. Madec)  add zps case 
     13   !!            3.0  ! 2010-06  (C. Ethe, G. Madec) Merge TRA-TRC 
     14   !!            3.7  ! 2014-01  (G. Madec, S. Masson) re-entrant laplacian  
    1415   !!---------------------------------------------------------------------- 
    1516 
    1617   !!---------------------------------------------------------------------- 
    17    !!   tra_ldf_lap  : update the tracer trend with the horizontal diffusion 
    18    !!                 using a iso-level harmonic (laplacien) operator. 
     18   !!   tra_ldf_lap : update the tracer trend with the lateral diffusion : iso-level laplacian operator 
     19   !!   tra_ldf_blp : update the tracer trend with the lateral diffusion : iso-level bilaplacian operator 
    1920   !!---------------------------------------------------------------------- 
    2021   USE oce             ! ocean dynamics and active tracers 
    2122   USE dom_oce         ! ocean space and time domain 
    22    USE ldftra_oce      ! ocean active tracers: lateral physics 
    23    USE in_out_manager  ! I/O manager 
     23   USE ldftra          ! lateral physics: eddy diffusivity 
    2424   USE diaptr          ! poleward transport diagnostics 
    2525   USE trc_oce         ! share passive tracers/Ocean variables 
    26    USE lib_mpp         ! MPP library 
     26   USE zpshde          ! partial step: hor. derivative     (zps_hde routine) 
     27   ! 
     28   USE in_out_manager  ! I/O manager 
     29   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     30   USE lib_mpp         ! distribued memory computing library 
    2731   USE timing          ! Timing 
     32   USE wrk_nemo        ! Memory allocation 
    2833 
    2934   IMPLICIT NONE 
    3035   PRIVATE 
    3136 
    32    PUBLIC   tra_ldf_lap   ! routine called by step.F90 
     37   PUBLIC   tra_ldf_lap   ! routine called by traldf.F90 
    3338 
    3439   !! * Substitutions 
    3540#  include "domzgr_substitute.h90" 
    36 #  include "ldftra_substitute.h90" 
    3741#  include "vectopt_loop_substitute.h90" 
    3842   !!---------------------------------------------------------------------- 
    39    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     43   !! NEMO/OPA 3.7 , NEMO Consortium (2014) 
    4044   !! $Id$ 
    4145   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    4347CONTAINS 
    4448 
    45    SUBROUTINE tra_ldf_lap( kt, kit000, cdtype, pgu , pgv ,    & 
    46       &                                        pgui, pgvi,    & 
    47       &                                ptb, pta, kjpt )  
     49   SUBROUTINE tra_ldf_lap( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   & 
     50      &                                                    pgui, pgvi,   & 
     51      &                                        ptb , pta , kjpt, kpass )  
    4852      !!---------------------------------------------------------------------- 
    4953      !!                  ***  ROUTINE tra_ldf_lap  *** 
    … …  
    5559      !!      fields (forward time scheme). The horizontal diffusive trends of  
    5660      !!      the tracer is given by: 
    57       !!          difft = 1/(e1t*e2t*e3t) {  di-1[ aht e2u*e3u/e1u di(tb) ] 
    58       !!                                   + dj-1[ aht e1v*e3v/e2v dj(tb) ] } 
     61      !!          difft = 1/(e1e2t*e3t) {  di-1[ pahu e2u*e3u/e1u di(tb) ] 
     62      !!                                 + dj-1[ pahv e1v*e3v/e2v dj(tb) ] } 
    5963      !!      Add this trend to the general tracer trend pta : 
    6064      !!          pta = pta + difft 
    … …  
    6367      !!                harmonic mixing trend. 
    6468      !!---------------------------------------------------------------------- 
    65       USE oce, ONLY:   ztu => ua , ztv => va  ! (ua,va) used as workspace 
    66       ! 
    6769      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index 
    68       INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index 
     70      INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index 
    6971      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
    7072      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers 
     73      INTEGER                              , INTENT(in   ) ::   kpass      ! =1/2 first or second passage 
     74      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s] 
    7175      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels 
    72       REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top levels 
     76      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels 
    7377      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields 
    7478      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend  
    7579      ! 
    76       INTEGER  ::   ji, jj, jk, jn       ! dummy loop indices 
    77       INTEGER  ::   iku, ikv, ierr       ! local integers 
    78       REAL(wp) ::   zabe1, zabe2, zbtr   ! local scalars 
     80      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
     81      REAL(wp) ::   zsign            ! local scalars 
     82      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztu, ztv, zaheeu, zaheev 
    7983      !!---------------------------------------------------------------------- 
    8084      ! 
    81       IF( nn_timing == 1 ) CALL timing_start('tra_ldf_lap') 
     85      IF( nn_timing == 1 )   CALL timing_start('tra_ldf_lap') 
    8286      ! 
    83       IF( kt == kit000 )  THEN 
    84          IF(lwp) WRITE(numout,*) 
    85          IF(lwp) WRITE(numout,*) 'tra_ldf_lap : iso-level laplacian diffusion on ', cdtype 
    86          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 
     87      IF( kt == nit000 .AND. lwp )  THEN 
     88         WRITE(numout,*) 
     89         WRITE(numout,*) 'tra_ldf_lap : iso-level laplacian diffusion on ', cdtype, ', pass=', kpass 
     90         WRITE(numout,*) '~~~~~~~~~~~ ' 
    8791      ENDIF 
    88  
    89       !                                                          ! =========== ! 
    90       DO jn = 1, kjpt                                            ! tracer loop ! 
    91          !                                                       ! =========== !     
    92          DO jk = 1, jpkm1                                            ! slab loop 
    93             !                                            
    94             ! 1. First derivative (gradient) 
    95             ! ------------------- 
     92      ! 
     93      CALL wrk_alloc( jpi,jpj,jpk,   ztu, ztv, zaheeu, zaheev )  
     94      ! 
     95      !                                !==  Initialization of metric arrays used for all tracers  ==! 
     96      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0) 
     97      ELSE                    ;   zsign = -1._wp 
     98      ENDIF 
     99      DO jk = 1, jpkm1 
     100         DO jj = 1, jpjm1 
     101            DO ji = 1, fs_jpim1   ! vector opt. 
     102               zaheeu(ji,jj,jk) = zsign * pahu(ji,jj,jk) * e2_e1u(ji,jj) * fse3u_n(ji,jj,jk)   !!gm   * umask(ji,jj,jk) pah masked! 
     103               zaheev(ji,jj,jk) = zsign * pahv(ji,jj,jk) * e1_e2v(ji,jj) * fse3v_n(ji,jj,jk)   !!gm   * vmask(ji,jj,jk) 
     104            END DO 
     105         END DO 
     106      END DO 
     107      ! 
     108      !                             ! =========== ! 
     109      DO jn = 1, kjpt               ! tracer loop ! 
     110         !                          ! =========== !     
     111         !                                
     112         DO jk = 1, jpkm1              !== First derivative (gradient)  ==! 
    96113            DO jj = 1, jpjm1 
    97                DO ji = 1, fs_jpim1   ! vector opt. 
    98                   zabe1 = fsahtu(ji,jj,jk) * umask(ji,jj,jk) * e2_e1u(ji,jj) * fse3u_n(ji,jj,jk) 
    99                   zabe2 = fsahtv(ji,jj,jk) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * fse3v_n(ji,jj,jk) 
    100                   ztu(ji,jj,jk) = zabe1 * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) 
    101                   ztv(ji,jj,jk) = zabe2 * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) 
     114               DO ji = 1, fs_jpim1 
     115                  ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) 
     116                  ztv(ji,jj,jk) = zaheev(ji,jj,jk) * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) 
    102117               END DO 
    103118            END DO 
    104             IF( ln_zps ) THEN      ! set gradient at partial step level for the last ocean cell 
     119         END DO   
     120         IF( ln_zps ) THEN                ! set gradient at bottom/top ocean level 
     121            DO jj = 1, jpjm1                    ! bottom 
     122               DO ji = 1, fs_jpim1 
     123                  ztu(ji,jj,mbku(ji,jj)) = zaheeu(ji,jj,mbku(ji,jj)) * pgu(ji,jj,jn) 
     124                  ztv(ji,jj,mbkv(ji,jj)) = zaheev(ji,jj,mbkv(ji,jj)) * pgv(ji,jj,jn) 
     125               END DO 
     126            END DO   
     127            IF( ln_isfcav ) THEN                ! top in ocean cavities only 
    105128               DO jj = 1, jpjm1 
    106129                  DO ji = 1, fs_jpim1   ! vector opt. 
    107                      ! last level 
    108                      iku = mbku(ji,jj) 
    109                      ikv = mbkv(ji,jj) 
    110                      IF( iku == jk ) THEN 
    111                         zabe1 = fsahtu(ji,jj,iku) * umask(ji,jj,iku) * e2_e1u(ji,jj) * fse3u_n(ji,jj,iku) 
    112                         ztu(ji,jj,jk) = zabe1 * pgu(ji,jj,jn) 
    113                      ENDIF 
    114                      IF( ikv == jk ) THEN 
    115                         zabe2 = fsahtv(ji,jj,ikv) * vmask(ji,jj,ikv) * e1_e2v(ji,jj) * fse3v_n(ji,jj,ikv) 
    116                         ztv(ji,jj,jk) = zabe2 * pgv(ji,jj,jn) 
    117                      ENDIF 
     130                     IF( miku(ji,jj) > 1 )   ztu(ji,jj,miku(ji,jj)) = zaheeu(ji,jj,miku(ji,jj)) * pgui(ji,jj,jn)  
     131                     IF( mikv(ji,jj) > 1 )   ztv(ji,jj,mikv(ji,jj)) = zaheev(ji,jj,mikv(ji,jj)) * pgvi(ji,jj,jn)  
    118132                  END DO 
    119133               END DO 
    120134            ENDIF 
    121             ! (ISH) 
    122             IF( ln_zps .AND. ln_isfcav ) THEN      ! set gradient at partial step level for the first ocean cell 
    123                                                    ! into a cavity 
    124                DO jj = 1, jpjm1 
    125                   DO ji = 1, fs_jpim1   ! vector opt. 
    126                      ! ice shelf level level MAX(2,jk) => only where ice shelf 
    127                      iku = miku(ji,jj)  
    128                      ikv = mikv(ji,jj)  
    129                      IF( iku == MAX(2,jk) ) THEN  
    130                         zabe1 = fsahtu(ji,jj,iku) * umask(ji,jj,iku) * e2_e1u(ji,jj) * fse3u_n(ji,jj,iku)  
    131                         ztu(ji,jj,jk) = zabe1 * pgui(ji,jj,jn)  
    132                      ENDIF  
    133                      IF( ikv == MAX(2,jk) ) THEN  
    134                         zabe2 = fsahtv(ji,jj,ikv) * vmask(ji,jj,ikv) * e1_e2v(ji,jj) * fse3v_n(ji,jj,ikv)  
    135                         ztv(ji,jj,jk) = zabe2 * pgvi(ji,jj,jn)  
    136                      END IF  
    137                   END DO 
    138                END DO 
    139             ENDIF 
    140           
    141           
    142             ! 2. Second derivative (divergence) added to the general tracer trends 
    143             ! --------------------------------------------------------------------- 
     135         ENDIF 
     136         ! 
     137         DO jk = 1, jpkm1              !== Second derivative (divergence) added to the general tracer trends  ==! 
    144138            DO jj = 2, jpjm1 
    145                DO ji = fs_2, fs_jpim1   ! vector opt. 
    146                   zbtr = 1._wp / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 
    147                   ! horizontal diffusive trends added to the general tracer trends 
    148                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   & 
    149                      &                                          + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) 
     139               DO ji = fs_2, fs_jpim1 
     140                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)     & 
     141                     &                                   + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) )   & 
     142                     &                                / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 
    150143               END DO 
    151144            END DO 
    152             ! 
    153          END DO                                             !  End of slab   
     145         END DO   
    154146         ! 
    155          ! "Poleward" diffusive heat or salt transports 
    156          IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
    157             IF( jn  == jp_tem)   htr_ldf(:) = ptr_sj( ztv(:,:,:) ) 
    158             IF( jn  == jp_sal)   str_ldf(:) = ptr_sj( ztv(:,:,:) ) 
     147         !                             !== "Poleward" diffusive heat or salt transports  ==! 
     148         IF( ( kpass == 1 .AND. .NOT.ln_traldf_blp ) .OR.  &     !==  first pass only (  laplacian)  ==! 
     149             ( kpass == 2 .AND.      ln_traldf_blp ) ) THEN      !==  2nd   pass only (bilaplacian)  ==! 
     150            IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
     151               IF( jn  == jp_tem)   htr_ldf(:) = ptr_sj( -ztv(:,:,:) ) 
     152               IF( jn  == jp_sal)   str_ldf(:) = ptr_sj( -ztv(:,:,:) ) 
     153            ENDIF 
    159154         ENDIF 
    160          !                                                  ! ================== 
    161       END DO                                                ! end of tracer loop 
    162       !                                                     ! ================== 
    163       IF( nn_timing == 1 ) CALL timing_stop('tra_ldf_lap') 
     155         !                          ! ================== 
     156      END DO                        ! end of tracer loop 
     157      !                             ! ================== 
     158      ! 
     159      CALL wrk_dealloc( jpi,jpj,jpk,   ztu, ztv, zaheeu, zaheev )  
     160      ! 
     161      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_lap') 
    164162      ! 
    165163   END SUBROUTINE tra_ldf_lap 
    166  
     164    
    167165   !!============================================================================== 
    168166END MODULE traldf_lap 

  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_triad.F90

    r5722 r5758  
    1 MODULE traldf_iso_grif 
     1MODULE traldf_triad 
    22   !!====================================================================== 
    3    !!                   ***  MODULE  traldf_iso_grif  *** 
     3   !!                   ***  MODULE  traldf_triad  *** 
    44   !! Ocean  tracers:  horizontal component of the lateral tracer mixing trend 
    55   !!====================================================================== 
    6    !! History : 3.3  ! 2010-10  (G. Nurser, C. Harris, G. Madec) 
    7    !!                !          Griffies operator version adapted from traldf_iso.F90 
     6   !! History :  3.3  ! 2010-10  (G. Nurser, C. Harris, G. Madec)  Griffies operator (original code) 
     7   !!            3.7  ! 2013-12  (F. Lemarie, G. Madec)  triad operator (Griffies) + Method of Stabilizing Correction 
    88   !!---------------------------------------------------------------------- 
    9 #if   defined key_ldfslp   ||   defined key_esopa 
     9 
    1010   !!---------------------------------------------------------------------- 
    11    !!   'key_ldfslp'               slope of the lateral diffusive direction 
    12    !!---------------------------------------------------------------------- 
    13    !!   tra_ldf_iso_grif  : update the tracer trend with the horizontal component 
    14    !!                       of the Griffies iso-neutral laplacian operator 
     11   !!   tra_ldf_triad : update the tracer trend with the iso-neutral laplacian triad-operator 
    1512   !!---------------------------------------------------------------------- 
    1613   USE oce             ! ocean dynamics and active tracers 
    … …  
    1916   USE trc_oce         ! share passive tracers/Ocean variables 
    2017   USE zdf_oce         ! ocean vertical physics 
    21    USE ldftra_oce      ! ocean active tracers: lateral physics 
    22    USE ldfslp          ! iso-neutral slopes 
     18   USE ldftra          ! lateral physics: eddy diffusivity 
     19   USE ldfslp          ! lateral physics: iso-neutral slopes 
     20   USE traldf_iso      ! lateral diffusion (Madec operator)         (tra_ldf_iso routine) 
    2321   USE diaptr          ! poleward transport diagnostics 
     22   USE zpshde          ! partial step: hor. derivative     (zps_hde routine) 
     23   ! 
    2424   USE in_out_manager  ! I/O manager 
    2525   USE iom             ! I/O library 
    … …  
    2929   USE timing          ! Timing 
    3030 
    31  
    3231   IMPLICIT NONE 
    3332   PRIVATE 
    3433 
    35    PUBLIC   tra_ldf_iso_grif   ! routine called by traldf.F90 
    36  
    37    REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE, SAVE ::   psix_eiv, psiy_eiv   !: eiv stream function (diag only) 
    38    REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE, SAVE ::   ah_wslp2             !: aeiv*w-slope^2 
    39    REAL(wp),         DIMENSION(:,:,:), ALLOCATABLE, SAVE ::   zdkt3d               !: vertical tracer gradient at 2 levels 
     34   PUBLIC   tra_ldf_triad   ! routine called by traldf.F90 
     35 
     36   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE ::   zdkt3d   !: vertical tracer gradient at 2 levels 
    4037 
    4138   !! * Substitutions 
    4239#  include "domzgr_substitute.h90" 
    43 #  include "ldftra_substitute.h90" 
    4440#  include "vectopt_loop_substitute.h90" 
    45 #  include "ldfeiv_substitute.h90" 
    4641   !!---------------------------------------------------------------------- 
    47    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     42   !! NEMO/OPA 3.7 , NEMO Consortium (2015) 
    4843   !! $Id$ 
    4944   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    5146CONTAINS 
    5247 
    53   SUBROUTINE tra_ldf_iso_grif( kt, kit000, cdtype, pgu, pgv,              & 
    54        &                                   ptb, pta, kjpt, pahtb0 ) 
     48  SUBROUTINE tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   & 
     49      &                                                     pgui, pgvi,   & 
     50      &                                         ptb , ptbb, pta , kjpt, kpass ) 
    5551      !!---------------------------------------------------------------------- 
    56       !!                  ***  ROUTINE tra_ldf_iso_grif  *** 
     52      !!                  ***  ROUTINE tra_ldf_triad  *** 
    5753      !! 
    5854      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive 
    … …  
    6662      !!      nal or geopotential slopes computed in routine ldfslp. 
    6763      !! 
    68       !!      1st part :  masked horizontal derivative of T  ( di[ t ] ) 
    69       !!      ========    with partial cell update if ln_zps=T. 
     64      !!      see documentation for the desciption 
    7065      !! 
    71       !!      2nd part :  horizontal fluxes of the lateral mixing operator 
    72       !!      ======== 
    73       !!         zftu = (aht+ahtb0) e2u*e3u/e1u di[ tb ] 
    74       !!               - aht       e2u*uslp    dk[ mi(mk(tb)) ] 
    75       !!         zftv = (aht+ahtb0) e1v*e3v/e2v dj[ tb ] 
    76       !!               - aht       e2u*vslp    dk[ mj(mk(tb)) ] 
    77       !!      take the horizontal divergence of the fluxes: 
    78       !!         difft = 1/(e1t*e2t*e3t) {  di-1[ zftu ] +  dj-1[ zftv ]  } 
    79       !!      Add this trend to the general trend (ta,sa): 
    80       !!         ta = ta + difft 
    81       !! 
    82       !!      3rd part: vertical trends of the lateral mixing operator 
    83       !!      ========  (excluding the vertical flux proportional to dk[t] ) 
    84       !!      vertical fluxes associated with the rotated lateral mixing: 
    85       !!         zftw =-aht {  e2t*wslpi di[ mi(mk(tb)) ] 
    86       !!                     + e1t*wslpj dj[ mj(mk(tb)) ]  } 
    87       !!      take the horizontal divergence of the fluxes: 
    88       !!         difft = 1/(e1t*e2t*e3t) dk[ zftw ] 
    89       !!      Add this trend to the general trend (ta,sa): 
    90       !!         pta = pta + difft 
    91       !! 
    92       !! ** Action :   Update pta arrays with the before rotated diffusion 
     66      !! ** Action :   pta   updated with the before rotated diffusion 
     67      !!               ah_wslp2 .... 
     68      !!               akz   stabilizing vertical diffusivity coefficient (used in trazdf_imp) 
    9369      !!---------------------------------------------------------------------- 
    94       USE oce     , ONLY:   zftu => ua       , zftv => va            ! (ua,va) used as 3D workspace 
    95       ! 
    9670      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index 
    9771      INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index 
    9872      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
    9973      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers 
     74      INTEGER                              , INTENT(in   ) ::   kpass      ! =1/2 first or second passage 
     75      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s] 
    10076      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels 
    101       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields 
     77      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels 
     78      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! tracer (kpass=1) or laplacian of tracer (kpass=2) 
     79      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptbb       ! tracer (only used in kpass=2) 
    10280      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend 
    103       REAL(wp)                             , INTENT(in   ) ::   pahtb0     ! background diffusion coef 
    104       ! 
    105       INTEGER  ::  ji, jj, jk,jn   ! dummy loop indices 
    106       INTEGER  ::  ip,jp,kp        ! dummy loop indices 
    107       INTEGER  ::  ierr            ! temporary integer 
    108       REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3   ! local scalars 
    109       REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4   !   -      - 
    110       REAL(wp) ::  zcoef0, zbtr                  !   -      - 
     81      ! 
     82      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
     83      INTEGER  ::  ip,jp,kp         ! dummy loop indices 
     84      INTEGER  ::  ierr            ! local integer 
     85      REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3          ! local scalars 
     86      REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4          !   -      - 
     87      REAL(wp) ::  zcoef0, ze3w_2, zsign, z2dt, z1_2dt  !   -      - 
    11188      ! 
    11289      REAL(wp) ::   zslope_skew, zslope_iso, zslope2, zbu, zbv 
    113       REAL(wp) ::   ze1ur, zdxt, ze2vr, ze3wr, zdyt, zdzt 
     90      REAL(wp) ::   ze1ur, ze2vr, ze3wr, zdxt, zdyt, zdzt 
    11491      REAL(wp) ::   zah, zah_slp, zaei_slp 
    115       REAL(wp), POINTER, DIMENSION(:,:  ) :: z2d 
    116       REAL(wp), POINTER, DIMENSION(:,:,:) :: zdit, zdjt, ztfw  
    117       REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d       ! 3D workspace 
     92#if defined key_diaar5 
     93      REAL(wp) ::   zztmp              ! local scalar 
     94#endif 
     95      REAL(wp), POINTER, DIMENSION(:,:  ) :: z2d                                            ! 2D workspace 
     96      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdit, zdjt, zftu, zftv, ztfw, zpsi_uw, zpsi_vw   ! 3D     - 
    11897      !!---------------------------------------------------------------------- 
    11998      ! 
    120       IF( nn_timing == 1 )  CALL timing_start('tra_ldf_iso_grif') 
    121       ! 
    122       CALL wrk_alloc( jpi, jpj,      z2d )  
    123       CALL wrk_alloc( jpi, jpj, jpk, zdit, zdjt, ztfw  )  
    124       ! 
    125  
    126       IF( kt == kit000 .AND. .NOT.ALLOCATED(ah_wslp2) )  THEN 
     99      IF( nn_timing == 1 )  CALL timing_start('tra_ldf_triad') 
     100      ! 
     101      CALL wrk_alloc( jpi,jpj,       z2d )  
     102      CALL wrk_alloc( jpi,jpj,jpk,   zdit, zdjt, zftu, zftv, ztfw, zpsi_uw, zpsi_vw  )  
     103      ! 
     104      IF( .NOT.ALLOCATED(zdkt3d) )  THEN 
     105         ALLOCATE( zdkt3d(jpi,jpj,0:1) , STAT=ierr ) 
     106         IF( lk_mpp   )   CALL mpp_sum ( ierr ) 
     107         IF( ierr > 0 )   CALL ctl_stop('STOP', 'tra_ldf_triad: unable to allocate arrays') 
     108      ENDIF 
     109     ! 
     110      IF( kpass == 1 .AND. kt == kit000 )  THEN 
    127111         IF(lwp) WRITE(numout,*) 
    128          IF(lwp) WRITE(numout,*) 'tra_ldf_iso_grif : rotated laplacian diffusion operator on ', cdtype 
    129          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    130          ALLOCATE( ah_wslp2(jpi,jpj,jpk) , zdkt3d(jpi,jpj,0:1), STAT=ierr ) 
    131          IF( lk_mpp   )   CALL mpp_sum ( ierr ) 
    132          IF( ierr > 0 )   CALL ctl_stop('STOP', 'tra_ldf_iso_grif: unable to allocate arrays') 
    133          IF( ln_traldf_gdia ) THEN 
    134             IF (.NOT. ALLOCATED(psix_eiv))THEN 
    135                 ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr ) 
    136                 IF( lk_mpp   )   CALL mpp_sum ( ierr ) 
    137                 IF( ierr > 0 )   CALL ctl_stop('STOP', 'tra_ldf_iso_grif: unable to allocate diagnostics') 
    138             ENDIF 
     112         IF(lwp) WRITE(numout,*) 'tra_ldf_triad : rotated laplacian diffusion operator on ', cdtype 
     113         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~' 
     114      ENDIF 
     115      ! 
     116      !                                               ! set time step size (Euler/Leapfrog) 
     117      IF( neuler == 0 .AND. kt == kit000 ) THEN   ;   z2dt =     rdttra(1)      ! at nit000   (Euler) 
     118      ELSE                                        ;   z2dt = 2.* rdttra(1)      !             (Leapfrog) 
     119      ENDIF 
     120      z1_2dt = 1._wp / z2dt 
     121      ! 
     122      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0) 
     123      ELSE                    ;   zsign = -1._wp 
     124      ENDIF 
     125                   
     126      !!---------------------------------------------------------------------- 
     127      !!   0 - calculate  ah_wslp2, akz, and optionally zpsi_uw, zpsi_vw 
     128      !!---------------------------------------------------------------------- 
     129      ! 
     130      IF( kpass == 1 ) THEN         !==  first pass only  and whatever the tracer is  ==! 
     131         ! 
     132         akz     (:,:,:) = 0._wp       
     133         ah_wslp2(:,:,:) = 0._wp 
     134         IF( ln_ldfeiv_dia ) THEN 
     135            zpsi_uw(:,:,:) = 0._wp 
     136            zpsi_vw(:,:,:) = 0._wp 
    139137         ENDIF 
    140       ENDIF 
    141  
    142       !!---------------------------------------------------------------------- 
    143       !!   0 - calculate  ah_wslp2, psix_eiv, psiy_eiv 
    144       !!---------------------------------------------------------------------- 
    145  
    146       !!gm Future development: consider using Ah defined at T-points and attached to the 4 t-point triads 
    147  
    148       ah_wslp2(:,:,:) = 0._wp 
    149       IF( ln_traldf_gdia ) THEN 
    150          psix_eiv(:,:,:) = 0._wp 
    151          psiy_eiv(:,:,:) = 0._wp 
    152       ENDIF 
    153  
    154       DO ip = 0, 1 
    155          DO kp = 0, 1 
    156             DO jk = 1, jpkm1 
    157                DO jj = 1, jpjm1 
    158                   DO ji = 1, fs_jpim1 
    159                      ze1ur = 1._wp / e1u(ji,jj) 
    160                      ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp) 
    161                      zbu   = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 
    162                      zah   = fsahtu(ji,jj,jk)                                  ! fsaht(ji+ip,jj,jk) 
    163                      zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 
    164                      ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 
    165                      ! (do this by *adding* gradient of depth) 
    166                      zslope2 = zslope_skew + ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * ze1ur * umask(ji,jj,jk+kp) 
    167                      zslope2 = zslope2 *zslope2 
    168                      ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp)    & 
    169                         &                     + zah * ( zbu * ze3wr / ( e1t(ji+ip,jj) * e2t(ji+ip,jj) ) ) * zslope2 
    170                      IF( ln_traldf_gdia ) THEN 
    171                         zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew           ! fsaeit(ji+ip,jj,jk)*zslope_skew 
    172                         psix_eiv(ji,jj,jk+kp) = psix_eiv(ji,jj,jk+kp) + 0.25_wp * zaei_slp 
    173                      ENDIF 
     138         ! 
     139         DO ip = 0, 1                            ! i-k triads 
     140            DO kp = 0, 1 
     141               DO jk = 1, jpkm1 
     142                  DO jj = 1, jpjm1 
     143                     DO ji = 1, fs_jpim1 
     144                        ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp) 
     145                        zbu   = e1e2u(ji,jj) * fse3u(ji,jj,jk) 
     146                        zah   = 0.25_wp * pahu(ji,jj,jk) 
     147                        zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 
     148                        ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth) 
     149                        zslope2 = zslope_skew + ( fsdept(ji+1,jj,jk) - fsdept(ji,jj,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp) 
     150                        zslope2 = zslope2 *zslope2 
     151                        ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji+ip,jj) * zslope2 
     152                        akz     (ji+ip,jj,jk+kp) = akz     (ji+ip,jj,jk+kp) + zah * r1_e1u(ji,jj)       & 
     153                           &                                                      * r1_e1u(ji,jj) * umask(ji,jj,jk+kp) 
     154                        ! 
     155                       IF( ln_ldfeiv_dia )   zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp)   & 
     156                           &                                       + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * zslope_skew 
     157                     END DO 
    174158                  END DO 
    175159               END DO 
    176160            END DO 
    177161         END DO 
    178       END DO 
    179       ! 
    180       DO jp = 0, 1 
    181          DO kp = 0, 1 
    182             DO jk = 1, jpkm1 
    183                DO jj = 1, jpjm1 
    184                   DO ji=1,fs_jpim1 
    185                      ze2vr = 1._wp / e2v(ji,jj) 
    186                      ze3wr = 1.0_wp / fse3w(ji,jj+jp,jk+kp) 
    187                      zbv   = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 
    188                      zah   = fsahtv(ji,jj,jk)                                  ! fsaht(ji,jj+jp,jk) 
    189                      zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 
    190                      ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 
    191                      !    (do this by *adding* gradient of depth) 
    192                      zslope2 = zslope_skew + ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) * ze2vr * vmask(ji,jj,jk+kp) 
    193                      zslope2 = zslope2 * zslope2 
    194                      ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp)   & 
    195                         &                     + zah * ( zbv * ze3wr / ( e1t(ji,jj+jp) * e2t(ji,jj+jp) ) ) * zslope2 
    196                      IF( ln_traldf_gdia ) THEN 
    197                         zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew           ! fsaeit(ji,jj+jp,jk)*zslope_skew 
    198                         psiy_eiv(ji,jj,jk+kp) = psiy_eiv(ji,jj,jk+kp) + 0.25_wp * zaei_slp 
    199                      ENDIF 
     162         ! 
     163         DO jp = 0, 1                            ! j-k triads  
     164            DO kp = 0, 1 
     165               DO jk = 1, jpkm1 
     166                  DO jj = 1, jpjm1 
     167                     DO ji = 1, fs_jpim1 
     168                        ze3wr = 1.0_wp / fse3w(ji,jj+jp,jk+kp) 
     169                        zbv   = e1e2v(ji,jj) * fse3v(ji,jj,jk) 
     170                        zah   = 0.25_wp * pahv(ji,jj,jk) 
     171                        zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 
     172                        ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 
     173                        !    (do this by *adding* gradient of depth) 
     174                        zslope2 = zslope_skew + ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp) 
     175                        zslope2 = zslope2 * zslope2 
     176                        ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj+jp) * zslope2 
     177                        akz     (ji,jj+jp,jk+kp) = akz     (ji,jj+jp,jk+kp) + zah * r1_e2v(ji,jj)     & 
     178                           &                                                      * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp) 
     179                        ! 
     180                        IF( ln_ldfeiv_dia )   zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp)   & 
     181                           &                                       + 0.25 * aeiv(ji,jj,jk) * e1v(ji,jj) * zslope_skew 
     182                     END DO 
    200183                  END DO 
    201184               END DO 
    202185            END DO 
    203186         END DO 
    204       END DO 
    205       ! 
    206       IF( iom_use("uoce_eiv") .OR. iom_use("voce_eiv") .OR. iom_use("woce_eiv") )  THEN 
    207          ! 
    208          IF( ln_traldf_gdia .AND. cdtype == 'TRA' ) THEN 
    209             CALL wrk_alloc( jpi , jpj , jpk  , zw3d ) 
    210             DO jk=1,jpkm1 
    211                zw3d(:,:,jk) = (psix_eiv(:,:,jk+1) - psix_eiv(:,:,jk))/fse3u(:,:,jk)  ! u_eiv = -dpsix/dz 
    212             END DO 
    213             zw3d(:,:,jpk) = 0._wp 
    214             CALL iom_put( "uoce_eiv", zw3d )    ! i-eiv current 
    215  
    216             DO jk=1,jpk-1 
    217                zw3d(:,:,jk) = (psiy_eiv(:,:,jk+1) - psiy_eiv(:,:,jk))/fse3v(:,:,jk)  ! v_eiv = -dpsiy/dz 
    218             END DO 
    219             zw3d(:,:,jpk) = 0._wp 
    220             CALL iom_put( "voce_eiv", zw3d )    ! j-eiv current 
    221  
    222             DO jk=1,jpk-1 
    223                DO jj = 2, jpjm1 
    224                   DO ji = fs_2, fs_jpim1  ! vector opt. 
    225                      zw3d(ji,jj,jk) = (psiy_eiv(ji,jj,jk) - psiy_eiv(ji,jj-1,jk))/e2t(ji,jj) + & 
    226                           &    (psix_eiv(ji,jj,jk) - psix_eiv(ji-1,jj,jk))/e1t(ji,jj) ! w_eiv = dpsiy/dy + dpsiy/dx 
    227                   END DO 
    228                END DO 
    229             END DO 
    230             zw3d(:,:,jpk) = 0._wp 
    231             CALL iom_put( "woce_eiv", zw3d )    ! vert. eiv current 
    232             CALL wrk_dealloc( jpi , jpj , jpk  , zw3d ) 
     187         ! 
     188         IF( ln_traldf_msc ) THEN                ! stabilizing vertical diffusivity coefficient 
     189            ! 
     190            IF( ln_traldf_blp ) THEN                ! bilaplacian operator 
     191               DO jk = 2, jpkm1 
     192                  DO jj = 1, jpjm1 
     193                     DO ji = 1, fs_jpim1 
     194                        akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   & 
     195                           &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( fse3w(ji,jj,jk) * fse3w(ji,jj,jk) )  ) 
     196                     END DO 
     197                  END DO 
     198               END DO 
     199            ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator 
     200               DO jk = 2, jpkm1 
     201                  DO jj = 1, jpjm1 
     202                     DO ji = 1, fs_jpim1 
     203                        ze3w_2 = fse3w(ji,jj,jk) * fse3w(ji,jj,jk) 
     204                        zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
     205                        akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt 
     206                     END DO 
     207                  END DO 
     208               END DO 
     209           ENDIF 
     210           ! 
     211         ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit 
     212            akz(:,:,:) = ah_wslp2(:,:,:)       
    233213         ENDIF 
    234214         ! 
    235       ENDIF 
    236       !                                                          ! =========== 
    237       DO jn = 1, kjpt                                            ! tracer loop 
    238          !                                                       ! =========== 
     215         IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' )   CALL ldf_eiv_dia( zpsi_uw, zpsi_vw ) 
     216         ! 
     217      ENDIF                                  !==  end 1st pass only  ==! 
     218      ! 
     219      !                                                           ! =========== 
     220      DO jn = 1, kjpt                                             ! tracer loop 
     221         !                                                        ! =========== 
    239222         ! Zero fluxes for each tracer 
     223!!gm  this should probably be done outside the jn loop 
    240224         ztfw(:,:,:) = 0._wp 
    241225         zftu(:,:,:) = 0._wp 
    242226         zftv(:,:,:) = 0._wp 
    243227         ! 
    244          DO jk = 1, jpkm1                          !==  before lateral T & S gradients at T-level jk  ==! 
     228         DO jk = 1, jpkm1        !==  before lateral T & S gradients at T-level jk  ==! 
    245229            DO jj = 1, jpjm1 
    246230               DO ji = 1, fs_jpim1   ! vector opt. 
    … …  
    250234            END DO 
    251235         END DO 
    252          IF( ln_zps.and.l_grad_zps ) THEN              ! partial steps: correction at the last level 
    253             DO jj = 1, jpjm1 
    254                DO ji = 1, jpim1 
     236         IF( ln_zps .AND. l_grad_zps ) THEN    ! partial steps: correction at top/bottom ocean level 
     237            DO jj = 1, jpjm1                       ! bottom level 
     238               DO ji = 1, fs_jpim1   ! vector opt. 
    255239                  zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 
    256240                  zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 
    257241               END DO 
    258242            END DO 
     243            IF( ln_isfcav ) THEN                   ! top level (ocean cavities only) 
     244               DO jj = 1, jpjm1 
     245                  DO ji = 1, fs_jpim1   ! vector opt. 
     246                     IF( miku(ji,jj)  > 1 )   zdit(ji,jj,miku(ji,jj) ) = pgui(ji,jj,jn)  
     247                     IF( mikv(ji,jj)  > 1 )   zdjt(ji,jj,mikv(ji,jj) ) = pgvi(ji,jj,jn)  
     248                  END DO 
     249               END DO 
     250            ENDIF 
    259251         ENDIF 
    260252 
    … …  
    272264            ELSE                 ;   zdkt3d(:,:,0) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * tmask(:,:,jk) 
    273265            ENDIF 
    274  
    275  
    276             IF (ln_botmix_grif) THEN 
     266            ! 
     267            zaei_slp = 0._wp 
     268            ! 
     269            IF( ln_botmix_triad ) THEN 
    277270               DO ip = 0, 1              !==  Horizontal & vertical fluxes 
    278271                  DO kp = 0, 1 
    279272                     DO jj = 1, jpjm1 
    280273                        DO ji = 1, fs_jpim1 
    281                            ze1ur = 1._wp / e1u(ji,jj) 
     274                           ze1ur = r1_e1u(ji,jj) 
     275                           zdxt  = zdit(ji,jj,jk) * ze1ur 
     276                           ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp) 
     277                           zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr 
     278                           zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 
     279                           zslope_iso  = triadi  (ji+ip,jj,jk,1-ip,kp) 
     280 
     281                           zbu = 0.25_wp * e1e2u(ji,jj) * fse3u(ji,jj,jk) 
     282                           ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????   ahu is masked.... 
     283                           zah = pahu(ji,jj,jk) 
     284                           zah_slp  = zah * zslope_iso 
     285                           IF( ln_ldfeiv )   zaei_slp = aeiu(ji,jj,jk) * zslope_skew 
     286                           zftu(ji   ,jj,jk   ) = zftu(ji   ,jj,jk   ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 
     287                           ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - ( zah_slp + zaei_slp) * zdxt                 * zbu * ze3wr 
     288                        END DO 
     289                     END DO 
     290                  END DO 
     291               END DO 
     292 
     293               DO jp = 0, 1 
     294                  DO kp = 0, 1 
     295                     DO jj = 1, jpjm1 
     296                        DO ji = 1, fs_jpim1 
     297                           ze2vr = r1_e2v(ji,jj) 
     298                           zdyt  = zdjt(ji,jj,jk) * ze2vr 
     299                           ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp) 
     300                           zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr 
     301                           zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 
     302                           zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp) 
     303                           zbv = 0.25_wp * e1e2v(ji,jj) * fse3v(ji,jj,jk) 
     304                           ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????  ahv is masked... 
     305                           zah = pahv(ji,jj,jk) 
     306                           zah_slp = zah * zslope_iso 
     307                           IF( ln_ldfeiv )   zaei_slp = aeiv(ji,jj,jk) * zslope_skew 
     308                           zftv(ji,jj   ,jk   ) = zftv(ji,jj   ,jk   ) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 
     309                           ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - ( zah_slp + zaei_slp ) * zdyt                * zbv * ze3wr 
     310                        END DO 
     311                     END DO 
     312                  END DO 
     313               END DO 
     314 
     315            ELSE 
     316 
     317               DO ip = 0, 1               !==  Horizontal & vertical fluxes 
     318                  DO kp = 0, 1 
     319                     DO jj = 1, jpjm1 
     320                        DO ji = 1, fs_jpim1 
     321                           ze1ur = r1_e1u(ji,jj) 
    282322                           zdxt  = zdit(ji,jj,jk) * ze1ur 
    283323                           ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp) 
    … …  
    286326                           zslope_iso  = triadi(ji+ip,jj,jk,1-ip,kp) 
    287327 
    288                            zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 
    289                            ! ln_botmix_grif is .T. don't mask zah for bottom half cells 
    290                            zah = fsahtu(ji,jj,jk)   !*umask(ji,jj,jk+kp)         !fsaht(ji+ip,jj,jk)           ===>>  ???? 
     328                           zbu = 0.25_wp * e1e2u(ji,jj) * fse3u(ji,jj,jk) 
     329                           ! ln_botmix_triad is .F. mask zah for bottom half cells 
     330                           zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp)         ! pahu(ji+ip,jj,jk)   ===>>  ???? 
    291331                           zah_slp  = zah * zslope_iso 
    292                            zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew    !fsaeit(ji+ip,jj,jk)*zslope_skew 
    293                            zftu(ji,jj,jk) = zftu(ji,jj,jk) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 
     332                           IF( ln_ldfeiv )   zaei_slp = aeiu(ji,jj,jk) * zslope_skew        ! fsaeit(ji+ip,jj,jk)*zslope_skew 
     333                           zftu(ji   ,jj,jk   ) = zftu(ji   ,jj,jk   ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 
    294334                           ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr 
    295335                        END DO 
    … …  
    302342                     DO jj = 1, jpjm1 
    303343                        DO ji = 1, fs_jpim1 
    304                            ze2vr = 1._wp / e2v(ji,jj) 
     344                           ze2vr = r1_e2v(ji,jj) 
    305345                           zdyt  = zdjt(ji,jj,jk) * ze2vr 
    306346                           ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp) 
    … …  
    308348                           zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 
    309349                           zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp) 
    310                            zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 
    311                            ! ln_botmix_grif is .T. don't mask zah for bottom half cells 
    312                            zah = fsahtv(ji,jj,jk)        !*vmask(ji,jj,jk+kp)  ! fsaht(ji,jj+jp,jk) 
     350                           zbv = 0.25_wp * e1e2v(ji,jj) * fse3v(ji,jj,jk) 
     351                           ! ln_botmix_triad is .F. mask zah for bottom half cells 
     352                           zah = pahv(ji,jj,jk) * vmask(ji,jj,jk+kp)         ! pahv(ji,jj+jp,jk)  ???? 
    313353                           zah_slp = zah * zslope_iso 
    314                            zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew        ! fsaeit(ji,jj+jp,jk)*zslope_skew 
     354                           IF( ln_ldfeiv )   zaei_slp = aeiv(ji,jj,jk) * zslope_skew        ! fsaeit(ji,jj+jp,jk)*zslope_skew 
    315355                           zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 
    316356                           ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr 
    … …  
    319359                  END DO 
    320360               END DO 
    321             ELSE 
    322                DO ip = 0, 1              !==  Horizontal & vertical fluxes 
    323                   DO kp = 0, 1 
    324                      DO jj = 1, jpjm1 
    325                         DO ji = 1, fs_jpim1 
    326                            ze1ur = 1._wp / e1u(ji,jj) 
    327                            zdxt  = zdit(ji,jj,jk) * ze1ur 
    328                            ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp) 
    329                            zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr 
    330                            zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 
    331                            zslope_iso  = triadi(ji+ip,jj,jk,1-ip,kp) 
    332  
    333                            zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 
    334                            ! ln_botmix_grif is .F. mask zah for bottom half cells 
    335                            zah = fsahtu(ji,jj,jk) * umask(ji,jj,jk+kp)         ! fsaht(ji+ip,jj,jk)   ===>>  ???? 
    336                            zah_slp  = zah * zslope_iso 
    337                            zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew        ! fsaeit(ji+ip,jj,jk)*zslope_skew 
    338                            zftu(ji,jj,jk) = zftu(ji,jj,jk) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 
    339                            ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr 
    340                         END DO 
    341                      END DO 
    342                   END DO 
    343                END DO 
    344  
    345                DO jp = 0, 1 
    346                   DO kp = 0, 1 
    347                      DO jj = 1, jpjm1 
    348                         DO ji = 1, fs_jpim1 
    349                            ze2vr = 1._wp / e2v(ji,jj) 
    350                            zdyt  = zdjt(ji,jj,jk) * ze2vr 
    351                            ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp) 
    352                            zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr 
    353                            zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 
    354                            zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp) 
    355                            zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 
    356                            ! ln_botmix_grif is .F. mask zah for bottom half cells 
    357                            zah = fsahtv(ji,jj,jk) * vmask(ji,jj,jk+kp)         ! fsaht(ji,jj+jp,jk) 
    358                            zah_slp = zah * zslope_iso 
    359                            zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew        ! fsaeit(ji,jj+jp,jk)*zslope_skew 
    360                            zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 
    361                            ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr 
    362                         END DO 
    363                      END DO 
    364                   END DO 
    365                END DO 
    366             END IF 
    367             !                          !==  divergence and add to the general trend  ==! 
     361            ENDIF 
     362            !                             !==  horizontal divergence and add to the general trend  ==! 
    368363            DO jj = 2 , jpjm1 
    369364               DO ji = fs_2, fs_jpim1  ! vector opt. 
    370                   zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    371                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zbtr * (   zftu(ji-1,jj,jk) - zftu(ji,jj,jk)   & 
    372                      &                                           + zftv(ji,jj-1,jk) - zftv(ji,jj,jk)   ) 
     365                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji-1,jj,jk) - zftu(ji,jj,jk)       & 
     366                     &                                           + zftv(ji,jj-1,jk) - zftv(ji,jj,jk)   )   & 
     367                     &                                        / (  e1e2t(ji,jj) * fse3t(ji,jj,jk)  ) 
    373368               END DO 
    374369            END DO 
    … …  
    376371         END DO 
    377372         ! 
    378          DO jk = 1, jpkm1              !== Divergence of vertical fluxes added to the general tracer trend 
     373         !                                !==  add the vertical 33 flux  ==! 
     374         IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz 
     375            DO jk = 2, jpkm1        
     376               DO jj = 1, jpjm1 
     377                  DO ji = fs_2, fs_jpim1 
     378                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk)   & 
     379                        &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             & 
     380                        &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 
     381                  END DO 
     382               END DO 
     383            END DO 
     384         ELSE                                   ! bilaplacian  
     385            SELECT CASE( kpass ) 
     386            CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2 
     387               DO jk = 2, jpkm1  
     388                  DO jj = 1, jpjm1 
     389                     DO ji = fs_2, fs_jpim1 
     390                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk)             & 
     391                           &                            * ah_wslp2(ji,jj,jk) * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 
     392                     END DO 
     393                  END DO 
     394               END DO  
     395            CASE(  2  )                            ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb  and ptbb gradients, resp. 
     396               DO jk = 2, jpkm1  
     397                  DO jj = 1, jpjm1 
     398                     DO ji = fs_2, fs_jpim1 
     399                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk)                      & 
     400                           &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   & 
     401                           &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   ) 
     402                     END DO 
     403                  END DO 
     404               END DO 
     405            END SELECT  
     406         ENDIF 
     407         ! 
     408         DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pta  ==! 
    379409            DO jj = 2, jpjm1 
    380410               DO ji = fs_2, fs_jpim1  ! vector opt. 
    381                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   & 
    382                      &                                / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     411                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   & 
     412                     &                                        / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    383413               END DO 
    384414            END DO 
    385415         END DO 
    386416         ! 
    387          !                             ! "Poleward" diffusive heat or salt transports (T-S case only) 
    388          IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
    389             IF( jn == jp_tem)   htr_ldf(:) = ptr_sj( zftv(:,:,:) )        ! 3.3  names 
    390             IF( jn == jp_sal)   str_ldf(:) = ptr_sj( zftv(:,:,:) ) 
    391          ENDIF 
    392  
    393          IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN 
    394            ! 
    395            IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN 
    396                z2d(:,:) = 0._wp  
    397                DO jk = 1, jpkm1 
    398                   DO jj = 2, jpjm1 
    399                      DO ji = fs_2, fs_jpim1   ! vector opt. 
    400                         z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk)  
    401                      END DO 
    402                   END DO 
    403                END DO 
    404                z2d(:,:) = rau0_rcp * z2d(:,:)  
    405                CALL lbc_lnk( z2d, 'U', -1. ) 
    406                CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction 
     417         IF( ( kpass == 1 .AND. ln_traldf_lap ) .OR.  &     !==  first pass only (  laplacian)  ==! 
     418             ( kpass == 2 .AND. ln_traldf_blp ) ) THEN      !==  2nd   pass      (bilaplacian)  ==! 
     419            ! 
     420            !                          ! "Poleward" diffusive heat or salt transports (T-S case only) 
     421            IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
     422               IF( jn == jp_tem)   htr_ldf(:) = ptr_sj( zftv(:,:,:) )        ! 3.3  names 
     423               IF( jn == jp_sal)   str_ldf(:) = ptr_sj( zftv(:,:,:) ) 
     424            ENDIF 
     425            ! 
     426            IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN 
     427              ! 
     428              IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN 
     429                  z2d(:,:) = zftu(ji,jj,1)  
     430                  DO jk = 2, jpkm1 
     431                     DO jj = 2, jpjm1 
     432                        DO ji = fs_2, fs_jpim1   ! vector opt. 
     433                           z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk)  
     434                        END DO 
     435                     END DO 
     436                  END DO 
     437                  z2d(:,:) = rau0_rcp * z2d(:,:)  
     438                  CALL lbc_lnk( z2d, 'U', -1. ) 
     439                  CALL iom_put( "udiff_heattr", z2d )                  ! heat i-transport 
     440                  ! 
     441                  z2d(:,:) = zftv(ji,jj,1)  
     442                  DO jk = 2, jpkm1 
     443                     DO jj = 2, jpjm1 
     444                        DO ji = fs_2, fs_jpim1   ! vector opt. 
     445                           z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk)  
     446                        END DO 
     447                     END DO 
     448                  END DO 
     449                  z2d(:,:) = rau0_rcp * z2d(:,:)      
     450                  CALL lbc_lnk( z2d, 'V', -1. ) 
     451                  CALL iom_put( "vdiff_heattr", z2d )                  !  heat j-transport 
     452               ENDIF 
    407453               ! 
    408                z2d(:,:) = 0._wp  
    409                DO jk = 1, jpkm1 
    410                   DO jj = 2, jpjm1 
    411                      DO ji = fs_2, fs_jpim1   ! vector opt. 
    412                         z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk)  
    413                      END DO 
    414                   END DO 
    415                END DO 
    416                z2d(:,:) = rau0_rcp * z2d(:,:)      
    417                CALL lbc_lnk( z2d, 'V', -1. ) 
    418                CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in i-direction 
    419             END IF 
    420             ! 
    421          ENDIF 
    422          ! 
    423       END DO 
    424       ! 
    425       CALL wrk_dealloc( jpi, jpj,      z2d )  
    426       CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt, ztfw  )  
    427       ! 
    428       IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_iso_grif') 
    429       ! 
    430   END SUBROUTINE tra_ldf_iso_grif 
    431  
    432 #else 
    433    !!---------------------------------------------------------------------- 
    434    !!   default option :   Dummy code   NO rotation of the diffusive tensor 
    435    !!---------------------------------------------------------------------- 
    436    REAL, PUBLIC, DIMENSION(:,:,:), ALLOCATABLE, SAVE ::   psix_eiv, psiy_eiv   !: eiv stream function (diag only) 
    437 CONTAINS 
    438    SUBROUTINE tra_ldf_iso_grif( kt, kit000, cdtype, pgu, pgv,              & 
    439        &                                   ptb, pta, kjpt, pahtb0 ) 
    440       CHARACTER(len=3) ::   cdtype 
    441       INTEGER          ::   kit000     ! first time step index 
    442       REAL, DIMENSION(:,:,:) ::   pgu, pgv   ! tracer gradient at pstep levels 
    443       REAL, DIMENSION(:,:,:,:) ::   ptb, pta 
    444       WRITE(*,*) 'tra_ldf_iso_grif: You should not have seen this print! error?', kt, cdtype,    & 
    445          &                  pgu(1,1,1), pgv(1,1,1), ptb(1,1,1,1), pta(1,1,1,1), kjpt, pahtb0 
    446    END SUBROUTINE tra_ldf_iso_grif 
    447 #endif 
     454            ENDIF 
     455            ! 
     456         ENDIF                                                    !== end pass selection  ==! 
     457         ! 
     458         !                                                        ! =============== 
     459      END DO                                                      ! end tracer loop 
     460      !                                                           ! =============== 
     461      ! 
     462      CALL wrk_dealloc( jpi,jpj,       z2d )  
     463      CALL wrk_dealloc( jpi,jpj,jpk,   zdit, zdjt, zftu, zftv, ztfw, zpsi_uw, zpsi_vw  )  
     464      ! 
     465      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_triad') 
     466      ! 
     467   END SUBROUTINE tra_ldf_triad 
    448468 
    449469   !!============================================================================== 
    450 END MODULE traldf_iso_grif 
     470END MODULE traldf_triad 

  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90

    r5656 r5758  
    3737   USE traqsr          ! penetrative solar radiation (needed for nksr) 
    3838   USE phycst          ! physical constant 
    39    USE ldftra_oce      ! lateral physics on tracers 
     39   USE ldftra          ! lateral physics on tracers 
     40   USE ldfslp 
    4041   USE bdy_oce         ! BDY open boundary condition variables 
    4142   USE bdytra          ! open boundary condition (bdy_tra routine) 

  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90

    r5385 r5758  
    1919   USE sbc_oce         ! surface boundary condition: ocean 
    2020   USE dynspg_oce 
     21   ! 
     22   USE ldftra          ! lateral diffusion: eddy diffusivity 
     23   USE ldfslp          ! lateral diffusion: iso-neutral slope  
    2124   USE trazdf_exp      ! vertical diffusion: explicit (tra_zdf_exp     routine) 
    2225   USE trazdf_imp      ! vertical diffusion: implicit (tra_zdf_imp     routine) 
    23    USE ldftra_oce      ! ocean active tracers: lateral physics 
     26   ! 
    2427   USE trd_oce         ! trends: ocean variables 
    2528   USE trdtra          ! trends manager: tracers  
    … …  
    4548#  include "vectopt_loop_substitute.h90" 
    4649   !!---------------------------------------------------------------------- 
    47    !! NEMO/OPA 3.7 , NEMO Consortium (2014) 
     50   !! NEMO/OPA 3.7 , NEMO Consortium (2015) 
    4851   !! $Id$ 
    4952   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    8891         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    8992      END SELECT 
     93!!gm WHY here !   and I don't like that ! 
    9094      ! DRAKKAR SSS control { 
    9195      ! JMM avoid negative salinities near river outlet ! Ugly fix 
    9296      ! JMM : restore negative salinities to small salinities: 
    9397      WHERE ( tsa(:,:,:,jp_sal) < 0._wp )   tsa(:,:,:,jp_sal) = 0.1_wp 
     98!!gm 
    9499 
    95100      IF( l_trdtra )   THEN                      ! save the vertical diffusive trends for further diagnostics 
    … …  
    98103            ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal) - tsb(:,:,jk,jp_sal) ) / r2dtra(jk) ) - ztrds(:,:,jk) 
    99104         END DO 
     105!!gm this should be moved in trdtra.F90 and done on all trends 
    100106         CALL lbc_lnk( ztrdt, 'T', 1. ) 
    101107         CALL lbc_lnk( ztrds, 'T', 1. ) 
     108!!gm 
    102109         CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt ) 
    103110         CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds ) 
    … …  
    123130      !!      nzdf = 0   explicit (time-splitting) scheme (ln_zdfexp=T) 
    124131      !!           = 1   implicit (euler backward) scheme (ln_zdfexp=F) 
    125       !!      NB: rotation of lateral mixing operator or TKE or KPP scheme, 
    126       !!      the implicit scheme is required. 
     132      !!      NB: rotation of lateral mixing operator or TKE & GLS schemes, 
     133      !!          an implicit scheme is required. 
    127134      !!---------------------------------------------------------------------- 
    128135      USE zdftke 
    129136      USE zdfgls 
    130       USE zdfkpp 
    131137      !!---------------------------------------------------------------------- 
    132138 
    … …  
    137143 
    138144      ! Force implicit schemes 
    139       IF( lk_zdftke .OR. lk_zdfgls .OR. lk_zdfkpp )   nzdf = 1      ! TKE, GLS or KPP physics 
    140       IF( ln_traldf_iso                           )   nzdf = 1      ! iso-neutral lateral physics 
    141       IF( ln_traldf_hor .AND. ln_sco              )   nzdf = 1      ! horizontal lateral physics in s-coordinate 
     145      IF( lk_zdftke .OR. lk_zdfgls   )   nzdf = 1   ! TKE, or GLS physics 
     146      IF( ln_traldf_iso              )   nzdf = 1   ! iso-neutral lateral physics 
     147      IF( ln_traldf_hor .AND. ln_sco )   nzdf = 1   ! horizontal lateral physics in s-coordinate 
    142148      IF( ln_zdfexp .AND. nzdf == 1 )   CALL ctl_stop( 'tra_zdf : If using the rotation of lateral mixing operator',   & 
    143             &                         ' TKE or KPP scheme, the implicit scheme is required, set ln_zdfexp = .false.' ) 
     149            &                         ' GLS or TKE scheme, the implicit scheme is required, set ln_zdfexp = .false.' ) 
    144150 
    145151      ! Test: esopa 

  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90

    r5120 r5758  
    1919   
    2020   !!---------------------------------------------------------------------- 
    21    !!   tra_zdf_imp : Update the tracer trend with the diagonal vertical   
    22    !!                 part of the mixing tensor. 
    23    !!---------------------------------------------------------------------- 
    24    USE oce             ! ocean dynamics and tracers variables 
    25    USE dom_oce         ! ocean space and time domain variables  
    26    USE zdf_oce         ! ocean vertical physics variables 
    27    USE trc_oce         ! share passive tracers/ocean variables 
    28    USE domvvl          ! variable volume 
    29    USE ldftra_oce      ! ocean active tracers: lateral physics 
    30    USE ldftra          ! lateral mixing type 
    31    USE ldfslp          ! lateral physics: slope of diffusion 
    32    USE zdfddm          ! ocean vertical physics: double diffusion 
    33    USE traldf_iso_grif ! active tracers: Griffies operator 
    34    USE in_out_manager  ! I/O manager 
    35    USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    36    USE lib_mpp         ! MPP library 
    37    USE wrk_nemo        ! Memory Allocation 
    38    USE timing          ! Timing 
     21   !!   tra_zdf_imp : Update the tracer trend with the diagonal vertical part of the mixing tensor. 
     22   !!---------------------------------------------------------------------- 
     23   USE oce              ! ocean dynamics and tracers variables 
     24   USE dom_oce          ! ocean space and time domain variables  
     25   USE zdf_oce          ! ocean vertical physics variables 
     26   USE trc_oce          ! share passive tracers/ocean variables 
     27   USE domvvl           ! variable volume 
     28   USE ldftra           ! lateral mixing type 
     29   USE ldfslp           ! lateral physics: slope of diffusion 
     30   USE zdfddm           ! ocean vertical physics: double diffusion 
     31   USE traldf_iso_triad ! active tracers: Method of Stabilizing Correction 
     32   ! 
     33   USE in_out_manager   ! I/O manager 
     34   USE lbclnk           ! ocean lateral boundary conditions (or mpp link) 
     35   USE lib_mpp          ! MPP library 
     36   USE wrk_nemo         ! Memory Allocation 
     37   USE timing           ! Timing 
    3938 
    4039   IMPLICIT NONE 
    … …  
    4746   !! * Substitutions 
    4847#  include "domzgr_substitute.h90" 
    49 #  include "ldftra_substitute.h90" 
    5048#  include "zdfddm_substitute.h90" 
    5149#  include "vectopt_loop_substitute.h90" 
    5250   !!---------------------------------------------------------------------- 
    53    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     51   !! NEMO/OPA 3.7 , NEMO Consortium (2015) 
    5452   !! $Id$ 
    5553   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    120118            ELSE                                            ;   zwt(:,:,2:jpk) = fsavs(:,:,2:jpk) 
    121119            ENDIF 
    122             DO jj=1, jpj 
    123                DO ji=1, jpi 
    124                   zwt(ji,jj,1) = 0._wp 
    125                END DO 
    126             END DO 
    127 ! 
    128 #if defined key_ldfslp 
    129             ! isoneutral diffusion: add the contribution  
    130             IF( ln_traldf_grif    ) THEN     ! Griffies isoneutral diff 
    131                DO jk = 2, jpkm1 
    132                   DO jj = 2, jpjm1 
    133                      DO ji = fs_2, fs_jpim1   ! vector opt. 
    134                         zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk)        
     120            zwt(:,:,1) = 0._wp 
     121            ! 
     122            IF( l_ldfslp ) THEN            ! isoneutral diffusion: add the contribution  
     123               IF( ln_traldf_msc  ) THEN     ! MSC iso-neutral operator  
     124                  DO jk = 2, jpkm1 
     125                     DO jj = 2, jpjm1 
     126                        DO ji = fs_2, fs_jpim1   ! vector opt. 
     127                           zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk)   
     128                        END DO 
    135129                     END DO 
    136130                  END DO 
    137                END DO 
    138             ELSE IF( l_traldf_rot ) THEN     ! standard isoneutral diff 
    139                DO jk = 2, jpkm1 
    140                   DO jj = 2, jpjm1 
    141                      DO ji = fs_2, fs_jpim1   ! vector opt. 
    142                         zwt(ji,jj,jk) = zwt(ji,jj,jk) + fsahtw(ji,jj,jk)                       & 
    143                            &                          * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   & 
    144                            &                             + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  ) 
     131               ELSE                          ! standard or triad iso-neutral operator 
     132                  DO jk = 2, jpkm1 
     133                     DO jj = 2, jpjm1 
     134                        DO ji = fs_2, fs_jpim1   ! vector opt. 
     135                           zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk) 
     136                        END DO 
    145137                     END DO 
    146138                  END DO 
    147                END DO 
     139               ENDIF 
    148140            ENDIF 
    149 #endif 
     141            ! 
    150142            ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked) 
    151143            DO jk = 1, jpkm1 
    … …  
    202194               ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,1) 
    203195               ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t(ji,jj,1) 
    204                pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn)                     & 
    205                   &                      + p2dt(1) * ze3tn * pta(ji,jj,1,jn) 
     196               pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) + p2dt(1) * ze3tn * pta(ji,jj,1,jn) 
    206197            END DO 
    207198         END DO 

  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90

    r5120 r5758  
    9393      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv  ! hor. grad of prd at u- & v-pts (bottom) 
    9494      ! 
    95       INTEGER  ::   ji, jj, jn      ! Dummy loop indices 
    96       INTEGER  ::   iku, ikv, ikum1, ikvm1   ! partial step level (ocean bottom level) at u- and v-points 
    97       REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv  ! temporary scalars 
    98       REAL(wp), DIMENSION(jpi,jpj)      ::  zri, zrj, zhi, zhj   ! NB: 3rd dim=1 to use eos 
    99       REAL(wp), DIMENSION(jpi,jpj,kjpt) ::  zti, ztj             !  
     95      INTEGER  ::   ji, jj, jn                  ! Dummy loop indices 
     96      INTEGER  ::   iku, ikv, ikum1, ikvm1      ! partial step level (ocean bottom level) at u- and v-points 
     97      REAL(wp) ::   ze3wu, ze3wv, zmaxu, zmaxv  ! local scalars 
     98      REAL(wp), DIMENSION(jpi,jpj)      ::   zri, zrj, zhi, zhj   ! NB: 3rd dim=1 to use eos 
     99      REAL(wp), DIMENSION(jpi,jpj,kjpt) ::   zti, ztj             !  
    100100      !!---------------------------------------------------------------------- 
    101101      ! 
    102       IF( nn_timing == 1 )  CALL timing_start( 'zps_hde') 
    103       ! 
    104       pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ; 
    105       zti (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ; 
    106       zhi (:,:  )=0.0_wp ; zhj (:,:  )=0.0_wp ; 
     102      IF( nn_timing == 1 )   CALL timing_start( 'zps_hde') 
     103      ! 
     104      pgtu(:,:,:)=0._wp   ;   zti (:,:,:)=0._wp   ;   zhi (:,:  )=0._wp 
     105      pgtv(:,:,:)=0._wp   ;   ztj (:,:,:)=0._wp   ;   zhj (:,:  )=0._wp 
    107106      ! 
    108107      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==! 
    … …  
    149148         ! 
    150149      END DO 
    151  
    152       ! horizontal derivative of density anomalies (rd) 
    153       IF( PRESENT( prd ) ) THEN         ! depth of the partial step level 
    154          pgru(:,:)=0.0_wp   ; pgrv(:,:)=0.0_wp ;  
     150      !                 
     151      IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part) 
     152         pgru(:,:) = 0._wp 
     153         pgrv(:,:) = 0._wp                ! depth of the partial step level 
    155154         DO jj = 1, jpjm1 
    156155            DO ji = 1, jpim1 
    … …  
    167166            END DO 
    168167         END DO 
    169  
    170          ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 
    171          ! step and store it in  zri, zrj for each  case 
    172          CALL eos( zti, zhi, zri )   
    173          CALL eos( ztj, zhj, zrj ) 
    174  
    175          ! Gradient of density at the last level  
    176          DO jj = 1, jpjm1 
     168         ! 
     169         CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj  
     170         CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj  
     171         ! 
     172         DO jj = 1, jpjm1                 ! Gradient of density at the last level  
    177173            DO ji = 1, jpim1 
    178174               iku = mbku(ji,jj) 
    … …  
    192188      END IF 
    193189      ! 
    194       IF( nn_timing == 1 )  CALL timing_stop( 'zps_hde') 
     190      IF( nn_timing == 1 )   CALL timing_stop( 'zps_hde') 
    195191      ! 
    196192   END SUBROUTINE zps_hde 
    197    ! 
    198    SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv,   & 
    199       &                          prd, pgru, pgrv, pmru, pmrv, pgzu, pgzv, pge3ru, pge3rv,  & 
    200       &                   pgtui, pgtvi, pgrui, pgrvi, pmrui, pmrvi, pgzui, pgzvi, pge3rui, pge3rvi ) 
     193 
     194 
     195   SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu , pgtv , pgtui, pgtvi,                                   & 
     196      &                              prd, pgru , pgrv , pmru , pmrv , pgzu , pgzv , pge3ru , pge3rv ,   & 
     197      &                                   pgrui, pgrvi, pmrui, pmrvi, pgzui, pgzvi, pge3rui, pge3rvi ) 
    201198      !!---------------------------------------------------------------------- 
    202199      !!                     ***  ROUTINE zps_hde  *** 
    … …  
    245242      !!              - pge3ru, pge3rv, pge3rui, pge3rvi: horizontal gradient of rho weighted by local e3w at u- & v-points  
    246243      !!---------------------------------------------------------------------- 
    247       INTEGER                              , INTENT(in   )           ::  kt          ! ocean time-step index 
    248       INTEGER                              , INTENT(in   )           ::  kjpt        ! number of tracers 
    249       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta         ! 4D tracers fields 
    250       REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv  ! hor. grad. of ptra at u- & v-pts  
    251       REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtui, pgtvi  ! hor. grad. of stra at u- & v-pts (ISF) 
    252       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields 
    253       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv      ! hor. grad of prd at u- & v-pts (bottom) 
    254       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmru, pmrv      ! hor. sum  of prd at u- & v-pts (bottom) 
    255       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzu, pgzv      ! hor. grad of z   at u- & v-pts (bottom) 
    256       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3ru, pge3rv  ! hor. grad of prd weighted by local e3w at u- & v-pts (bottom) 
    257       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgrui, pgrvi      ! hor. grad of prd at u- & v-pts (top) 
    258       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmrui, pmrvi      ! hor. sum  of prd at u- & v-pts (top) 
    259       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzui, pgzvi      ! hor. grad of z   at u- & v-pts (top) 
    260       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3rui, pge3rvi  ! hor. grad of prd weighted by local e3w at u- & v-pts (top) 
     244      INTEGER                              , INTENT(in   )           ::  kt                ! ocean time-step index 
     245      INTEGER                              , INTENT(in   )           ::  kjpt              ! number of tracers 
     246      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta               ! 4D tracers fields 
     247      !                                                              !!  u-point ! v-point ! 
     248      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu    , pgtv    ! bottom GRADh( ptra )   
     249      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtui   , pgtvi   ! top    GRADh( ptra ) 
     250      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd               ! 3D density anomaly fields 
     251      !                                                              !!  u-point ! v-point ! 
     252      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru    , pgrv    ! bottom GRADh( prd  ) 
     253      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmru    , pmrv    ! bottom SUM  ( prd  ) 
     254      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzu    , pgzv    ! bottom GRADh( z    )  
     255      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3ru  , pge3rv  ! bottom GRADh( prd  ) weighted by e3w 
     256      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgrui   , pgrvi   ! top    GRADh( prd  )  
     257      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmrui   , pmrvi   ! top    SUM  ( prd  )  
     258      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzui   , pgzvi   ! top    GRADh( z    )  
     259      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3rui , pge3rvi ! top    GRADh( prd  ) weighted by e3w 
    261260      ! 
    262261      INTEGER  ::   ji, jj, jn      ! Dummy loop indices 
    … …  
    269268      IF( nn_timing == 1 )  CALL timing_start( 'zps_hde_isf') 
    270269      ! 
    271       pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ; 
    272       pgtui(:,:,:)=0.0_wp ; pgtvi(:,:,:)=0.0_wp ; 
    273       zti (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ; 
    274       zhi (:,:  )=0.0_wp ; zhj (:,:  )=0.0_wp ; 
     270      pgtu (:,:,:) = 0._wp   ;   pgtv (:,:,:) =0._wp 
     271      pgtui(:,:,:) = 0._wp   ;   pgtvi(:,:,:) =0._wp 
     272      zti  (:,:,:) = 0._wp   ;   ztj  (:,:,:) =0._wp 
     273      zhi  (:,:  ) = 0._wp   ;   zhj  (:,:  ) =0._wp 
    275274      ! 
    276275      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==! 
    … …  
    322321      END DO 
    323322 
    324       ! horizontal derivative of density anomalies (rd) 
    325       IF( PRESENT( prd ) ) THEN         ! depth of the partial step level 
    326          pgru(:,:)=0.0_wp   ; pgrv(:,:)=0.0_wp ;  
    327          pgzu(:,:)=0.0_wp   ; pgzv(:,:)=0.0_wp ; 
    328          pmru(:,:)=0.0_wp   ; pmru(:,:)=0.0_wp ; 
    329          pge3ru(:,:)=0.0_wp ; pge3rv(:,:)=0.0_wp ; 
    330          DO jj = 1, jpjm1 
     323      IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part) 
     324         ! 
     325         pgru  (:,:)=0._wp   ;   pgrv  (:,:) = 0._wp 
     326         pgzu  (:,:)=0._wp   ;   pgzv  (:,:) = 0._wp  
     327         pmru  (:,:)=0._wp   ;   pmru  (:,:) = 0._wp  
     328         pge3ru(:,:)=0._wp   ;   pge3rv(:,:) = 0._wp  
     329         ! 
     330         DO jj = 1, jpjm1                 ! depth of the partial step level 
    331331            DO ji = 1, jpim1 
    332332               iku = mbku(ji,jj) 
    … …  
    334334               ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 
    335335               ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 
    336  
     336               ! 
    337337               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) - ze3wu     ! i-direction: case 1 
    338338               ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) + ze3wu    ! -     -      case 2 
    … …  
    343343            END DO 
    344344         END DO 
    345           
    346          ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 
    347          ! step and store it in  zri, zrj for each  case 
    348          CALL eos( zti, zhi, zri )   
    349          CALL eos( ztj, zhj, zrj ) 
    350  
    351          ! Gradient of density at the last level  
    352          DO jj = 1, jpjm1 
     345         ! 
     346         CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj  
     347         CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj  
     348 
     349         DO jj = 1, jpjm1                 ! Gradient of density at the last level  
    353350            DO ji = 1, jpim1 
    354351               iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
    … …  
    394391         ! 
    395392      END IF 
    396          ! (ISH)  compute grui and gruvi 
     393      ! 
     394      !     !==  (ISH)  compute grui and gruvi  ==! 
     395      ! 
    397396      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!            ! 
    398397         DO jj = 1, jpjm1 
    … …  
    442441      END DO 
    443442 
    444       ! horizontal derivative of density anomalies (rd) 
    445       IF( PRESENT( prd ) ) THEN         ! depth of the partial step level 
     443      IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part) 
     444         ! 
    446445         pgrui(:,:)  =0.0_wp ; pgrvi(:,:)  =0.0_wp ; 
    447446         pgzui(:,:)  =0.0_wp ; pgzvi(:,:)  =0.0_wp ; 
    448447         pmrui(:,:)  =0.0_wp ; pmrui(:,:)  =0.0_wp ; 
    449448         pge3rui(:,:)=0.0_wp ; pge3rvi(:,:)=0.0_wp ; 
    450  
    451          DO jj = 1, jpjm1 
     449         ! 
     450         DO jj = 1, jpjm1        ! depth of the partial step level 
    452451            DO ji = 1, jpim1 
    453452               iku = miku(ji,jj) 
    … …  
    455454               ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 
    456455               ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 
    457  
     456               ! 
    458457               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) + ze3wu    ! i-direction: case 1 
    459458               ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) - ze3wu    ! -     -      case 2 
    … …  
    464463            END DO 
    465464         END DO 
    466  
    467          ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 
    468          ! step and store it in  zri, zrj for each  case 
    469          CALL eos( zti, zhi, zri )   
    470          CALL eos( ztj, zhj, zrj ) 
    471  
    472          ! Gradient of density at the last level  
    473          DO jj = 1, jpjm1 
     465         ! 
     466         CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj  
     467         CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj  
     468         ! 
     469         DO jj = 1, jpjm1                 ! Gradient of density at the last level  
    474470            DO ji = 1, jpim1 
    475471               iku = miku(ji,jj) ; ikup1 = miku(ji,jj) + 1 
    … …  
    482478                 pmrui  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) + prd(ji,jj,iku) )          ! i: 1  
    483479                 pge3rui(ji,jj) = umask(ji,jj,iku+1)                                                                  & 
    484                                 * ( (fse3w(ji+1,jj,iku+1) - ze3wu) * (zri(ji,jj    ) + prd(ji+1,jj,iku+1) + 2._wp)   & 
    485                                    - fse3w(ji  ,jj,iku+1)          * (prd(ji,jj,iku) + prd(ji  ,jj,iku+1) + 2._wp)   ) ! i: 1 
     480                    &           * ( (fse3w(ji+1,jj,iku+1) - ze3wu) * (zri(ji,jj    ) + prd(ji+1,jj,iku+1) + 2._wp)   & 
     481                    &              - fse3w(ji  ,jj,iku+1)          * (prd(ji,jj,iku) + prd(ji  ,jj,iku+1) + 2._wp)   ) ! i: 1 
    486482               ELSE 
    487483                 pgzui  (ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu) 
    … …  
    489485                 pmrui  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) + zri(ji,jj) )      ! i: 2 
    490486                 pge3rui(ji,jj) = umask(ji,jj,iku+1)                                                                   & 
    491                                 * (  fse3w(ji+1,jj,iku+1)          * (prd(ji+1,jj,iku) + prd(ji+1,jj,iku+1) + 2._wp)  & 
    492                                    -(fse3w(ji  ,jj,iku+1) + ze3wu) * (zri(ji,jj      ) + prd(ji  ,jj,iku+1) + 2._wp)  )     ! i: 2 
     487                    &           * (  fse3w(ji+1,jj,iku+1)          * (prd(ji+1,jj,iku) + prd(ji+1,jj,iku+1) + 2._wp)  & 
     488                    &              -(fse3w(ji  ,jj,iku+1) + ze3wu) * (zri(ji,jj      ) + prd(ji  ,jj,iku+1) + 2._wp)  )     ! i: 2 
    493489               ENDIF 
    494490               IF( ze3wv >= 0._wp ) THEN 
    … …  
    497493                 pmrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )        ! j: 1 
    498494                 pge3rvi(ji,jj) = vmask(ji,jj,ikv+1)                                                                  &  
    499                                 * ( (fse3w(ji,jj+1,ikv+1) - ze3wv) * ( zrj(ji,jj    ) + prd(ji,jj+1,ikv+1) + 2._wp)  & 
    500                                    - fse3w(ji,jj  ,ikv+1)          * ( prd(ji,jj,ikv) + prd(ji,jj  ,ikv+1) + 2._wp)  ) ! j: 1 
     495                     &           * ( (fse3w(ji,jj+1,ikv+1) - ze3wv) * ( zrj(ji,jj    ) + prd(ji,jj+1,ikv+1) + 2._wp)  & 
     496                     &                - fse3w(ji,jj  ,ikv+1)          * ( prd(ji,jj,ikv) + prd(ji,jj  ,ikv+1) + 2._wp)  ) ! j: 1 
    501497                                  ! + 2 due to the formulation in density and not in anomalie in hpg sco 
    502498               ELSE 
    … …  
    505501                 pmrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )     ! j: 2 
    506502                 pge3rvi(ji,jj) = vmask(ji,jj,ikv+1)                                                                   & 
    507                                 * (  fse3w(ji,jj+1,ikv+1)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikv+1) + 2._wp) & 
    508                                    -(fse3w(ji,jj  ,ikv+1) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikv+1) + 2._wp) )  ! j: 2 
     503                    &           * (  fse3w(ji,jj+1,ikv+1)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikv+1) + 2._wp) & 
     504                    &              -(fse3w(ji,jj  ,ikv+1) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikv+1) + 2._wp) )  ! j: 2 
    509505               ENDIF 
    510506            END DO 
    … …  
    517513      END IF   
    518514      ! 
    519       IF( nn_timing == 1 )  CALL timing_stop( 'zps_hde_isf') 
     515      IF( nn_timing == 1 )   CALL timing_stop( 'zps_hde_isf') 
    520516      ! 
    521517   END SUBROUTINE zps_hde_isf 
Note: See TracChangeset for help on using the changeset viewer.