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 5836 for trunk/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Ignore:
Timestamp:
2015-10-26T15:49:40+01:00 (9 years ago)
Author:
cetlod
Message:

merge the simplification branch onto the trunk, see ticket #1612

Location:
trunk/NEMOGCM/NEMO/OPA_SRC/TRA
Files:
8 deleted
14 edited
5 copied

Legend:

Unmodified
Added
Removed

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90

    r5147 r5836  
    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 
    9    !!---------------------------------------------------------------------- 
    10  
    11    !!---------------------------------------------------------------------- 
    12    !!   tra_adv      : compute ocean tracer advection trend 
    13    !!   tra_adv_ctl  : control the different options of advection scheme 
    14    !!---------------------------------------------------------------------- 
    15    USE oce             ! ocean dynamics and active tracers 
    16    USE dom_oce         ! ocean space and time domain 
    17    USE domvvl          ! variable vertical scale factors 
    18    USE traadv_cen2     ! 2nd order centered scheme (tra_adv_cen2   routine) 
    19    USE traadv_tvd      ! TVD      scheme           (tra_adv_tvd    routine) 
    20    USE traadv_muscl    ! MUSCL    scheme           (tra_adv_muscl  routine) 
    21    USE traadv_muscl2   ! MUSCL2   scheme           (tra_adv_muscl2 routine) 
    22    USE traadv_ubs      ! UBS      scheme           (tra_adv_ubs    routine) 
    23    USE traadv_qck      ! QUICKEST scheme           (tra_adv_qck    routine) 
    24    USE traadv_eiv      ! eddy induced velocity     (tra_adv_eiv    routine) 
    25    USE traadv_mle      ! ML eddy induced velocity  (tra_adv_mle    routine) 
    26    USE cla             ! cross land advection      (cla_traadv     routine) 
    27    USE ldftra_oce      ! lateral diffusion coefficient on tracers 
     8   !!            3.6  !  2011-06  (G. Madec)  Addition of Mixed Layer Eddy parameterisation 
     9   !!            3.7  !  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 
     11   !!---------------------------------------------------------------------- 
     12 
     13   !!---------------------------------------------------------------------- 
     14   !!   tra_adv       : compute ocean tracer advection trend 
     15   !!   tra_adv_ctl   : control the different options of advection scheme 
     16   !!---------------------------------------------------------------------- 
     17   USE oce            ! ocean dynamics and active tracers 
     18   USE dom_oce        ! ocean space and time domain 
     19   USE domvvl         ! variable vertical scale factors 
     20   USE traadv_cen     ! centered scheme           (tra_adv_cen  routine) 
     21   USE traadv_fct     ! FCT      scheme           (tra_adv_fct routine) 
     22   USE traadv_mus     ! MUSCL    scheme           (tra_adv_mus  routine) 
     23   USE traadv_ubs     ! UBS      scheme           (tra_adv_ubs  routine) 
     24   USE traadv_qck     ! QUICKEST scheme           (tra_adv_qck  routine) 
     25   USE traadv_mle     ! ML eddy induced velocity  (tra_adv_mle  routine) 
     26   USE ldftra         ! lateral diffusion: eddy diffusivity & EIV coeff. 
     27   USE ldfslp         ! Lateral diffusion: slopes of neutral surfaces 
    2828   ! 
    29    USE in_out_manager  ! I/O manager 
    30    USE iom             ! I/O module 
    31    USE prtctl          ! Print control 
    32    USE lib_mpp         ! MPP library 
    33    USE wrk_nemo        ! Memory Allocation 
    34    USE timing          ! Timing 
    35    USE sbc_oce 
     29   USE in_out_manager ! I/O manager 
     30   USE iom            ! I/O module 
     31   USE prtctl         ! Print control 
     32   USE lib_mpp        ! MPP library 
     33   USE wrk_nemo       ! Memory Allocation 
     34   USE timing         ! Timing 
     35 
    3636   USE diaptr          ! Poleward heat transport  
    37  
    3837 
    3938   IMPLICIT NONE 
    … …  
    4342   PUBLIC   tra_adv_init   ! routine called by opa module 
    4443 
    45    !                              !!* Namelist namtra_adv * 
    46    LOGICAL ::   ln_traadv_cen2     ! 2nd order centered scheme flag 
    47    LOGICAL ::   ln_traadv_tvd      ! TVD scheme flag 
    48    LOGICAL ::   ln_traadv_tvd_zts  ! TVD scheme flag with vertical sub time-stepping 
    49    LOGICAL ::   ln_traadv_muscl    ! MUSCL scheme flag 
    50    LOGICAL ::   ln_traadv_muscl2   ! MUSCL2 scheme flag 
    51    LOGICAL ::   ln_traadv_ubs      ! UBS scheme flag 
    52    LOGICAL ::   ln_traadv_qck      ! QUICKEST scheme flag 
    53    LOGICAL ::   ln_traadv_msc_ups  ! use upstream scheme within muscl 
    54  
    55  
    56    INTEGER ::   nadv   ! choice of the type of advection scheme 
    57  
     44   !                            !!* Namelist namtra_adv * 
     45   LOGICAL ::   ln_traadv_cen    ! centered scheme flag 
     46   INTEGER ::      nn_cen_h, nn_cen_v   ! =2/4 : horizontal and vertical choices of the order of CEN scheme 
     47   LOGICAL ::   ln_traadv_fct    ! FCT scheme flag 
     48   INTEGER ::      nn_fct_h, nn_fct_v   ! =2/4 : horizontal and vertical choices of the order of FCT scheme 
     49   INTEGER ::      nn_fct_zts           ! >=1 : 2nd order FCT with vertical sub-timestepping 
     50   LOGICAL ::   ln_traadv_mus    ! MUSCL scheme flag 
     51   LOGICAL ::      ln_mus_ups           ! use upstream scheme in vivcinity of river mouths 
     52   LOGICAL ::   ln_traadv_ubs    ! UBS scheme flag 
     53   INTEGER ::      nn_ubs_v             ! =2/4 : vertical choice of the order of UBS scheme 
     54   LOGICAL ::   ln_traadv_qck    ! QUICKEST scheme flag 
     55 
     56   INTEGER ::              nadv             ! choice of the type of advection scheme 
     57   ! 
     58   !                                        ! associated indices: 
     59   INTEGER, PARAMETER ::   np_NO_adv  = 0   ! no T-S advection 
     60   INTEGER, PARAMETER ::   np_CEN     = 1   ! 2nd/4th order centered scheme 
     61   INTEGER, PARAMETER ::   np_FCT     = 2   ! 2nd/4th order Flux Corrected Transport scheme 
     62   INTEGER, PARAMETER ::   np_FCT_zts = 3   ! 2nd order FCT scheme with vertical sub-timestepping 
     63   INTEGER, PARAMETER ::   np_MUS     = 4   ! MUSCL scheme 
     64   INTEGER, PARAMETER ::   np_UBS     = 5   ! 3rd order Upstream Biased Scheme 
     65   INTEGER, PARAMETER ::   np_QCK     = 6   ! QUICK scheme 
     66    
    5867   !! * Substitutions 
    5968#  include "domzgr_substitute.h90" 
    6069#  include "vectopt_loop_substitute.h90" 
    6170   !!---------------------------------------------------------------------- 
    62    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     71   !! NEMO/OPA 3.7 , NEMO Consortium (2014) 
    6372   !! $Id$ 
    6473   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    7483      !! ** Method  : - Update (ua,va) with the advection term following nadv 
    7584      !!---------------------------------------------------------------------- 
    76       ! 
    7785      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    7886      ! 
    … …  
    8391      IF( nn_timing == 1 )  CALL timing_start('tra_adv') 
    8492      ! 
    85       CALL wrk_alloc( jpi, jpj, jpk, zun, zvn, zwn ) 
     93      CALL wrk_alloc( jpi,jpj,jpk,   zun, zvn, zwn ) 
     94      ! 
    8695      !                                          ! set time step 
    8796      IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000 
    … …  
    91100      ENDIF 
    92101      ! 
    93       IF( nn_cla == 1 .AND. cp_cfg == 'orca' .AND. jp_cfg == 2 )   CALL cla_traadv( kt )       !==  Cross Land Advection  ==! (hor. advection) 
    94       ! 
    95       !                                               !==  effective transport  ==! 
     102      !                                         !==  effective transport  ==! 
    96103      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) 
     104         zun(:,:,jk) = e2u  (:,:) * fse3u(:,:,jk) * un(:,:,jk)                  ! eulerian transport only 
     105         zvn(:,:,jk) = e1v  (:,:) * fse3v(:,:,jk) * vn(:,:,jk) 
     106         zwn(:,:,jk) = e1e2t(:,:)                 * wn(:,:,jk) 
    100107      END DO 
    101108      ! 
    102       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 
     109      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN                                ! add z-tilde and/or vvl corrections 
    103110         zun(:,:,:) = zun(:,:,:) + un_td(:,:,:) 
    104111         zvn(:,:,:) = zvn(:,:,:) + vn_td(:,:,:) 
    105112      ENDIF 
    106113      ! 
    107       zun(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom 
    108       zvn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom 
    109       zwn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom 
    110       ! 
    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       
     114      zun(:,:,jpk) = 0._wp                                                      ! no transport trough the bottom 
     115      zvn(:,:,jpk) = 0._wp 
     116      zwn(:,:,jpk) = 0._wp 
     117      ! 
     118      IF( ln_ldfeiv .AND. .NOT. ln_traldf_triad )   & 
     119         &              CALL ldf_eiv_trp( kt, nit000, zun, zvn, zwn, 'TRA' )   ! add the eiv transport (if necessary) 
     120      ! 
     121      IF( ln_mle    )   CALL tra_adv_mle( kt, nit000, zun, zvn, zwn, 'TRA' )   ! add the mle transport (if necessary) 
     122      ! 
     123      CALL iom_put( "uocetr_eff", zun )                                        ! output effective transport       
    117124      CALL iom_put( "vocetr_eff", zvn ) 
    118125      CALL iom_put( "wocetr_eff", zwn ) 
    119126      ! 
    120       IF( ln_diaptr )   CALL dia_ptr( zvn )                                     ! diagnose the effective MSF  
    121       ! 
    122     
    123       SELECT CASE ( nadv )                            !==  compute advection trend and add it to general trend  ==! 
    124       CASE ( 1 )   ;    CALL tra_adv_cen2   ( kt, nit000, 'TRA',         zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  2nd order centered 
    125       CASE ( 2 )   ;    CALL tra_adv_tvd    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  TVD  
    126       CASE ( 3 )   ;    CALL tra_adv_muscl  ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb,      tsa, jpts, ln_traadv_msc_ups )   !  MUSCL  
    127       CASE ( 4 )   ;    CALL tra_adv_muscl2 ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  MUSCL2  
    128       CASE ( 5 )   ;    CALL tra_adv_ubs    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  UBS  
    129       CASE ( 6 )   ;    CALL tra_adv_qck    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  QUICKEST  
    130       CASE ( 7 )   ;    CALL tra_adv_tvd_zts( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  TVD ZTS 
    131       ! 
    132       CASE (-1 )                                      !==  esopa: test all possibility with control print  ==! 
    133          CALL tra_adv_cen2  ( kt, nit000, 'TRA',         zun, zvn, zwn, tsb, tsn, tsa, jpts )           
    134          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv0 - Ta: ', mask1=tmask,               & 
    135             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    136          CALL tra_adv_tvd   ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )           
    137          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv1 - Ta: ', mask1=tmask,               & 
    138             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    139          CALL tra_adv_muscl ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb,      tsa, jpts, ln_traadv_msc_ups )           
    140          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv3 - Ta: ', mask1=tmask,               & 
    141             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    142          CALL tra_adv_muscl2( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )           
    143          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv4 - Ta: ', mask1=tmask,               & 
    144             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    145          CALL tra_adv_ubs   ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )           
    146          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv5 - Ta: ', mask1=tmask,               & 
    147             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    148          CALL tra_adv_qck   ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )           
    149          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv6 - Ta: ', mask1=tmask,               & 
    150             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
     127!!gm ??? 
     128      IF( ln_diaptr )   CALL dia_ptr( zvn )                                    ! diagnose the effective MSF  
     129!!gm ??? 
     130      ! 
     131      SELECT CASE ( nadv )                      !==  compute advection trend and add it to general trend  ==! 
     132      ! 
     133      CASE ( np_CEN )                                    ! Centered scheme : 2nd / 4th order 
     134         CALL tra_adv_cen    ( kt, nit000, 'TRA',         zun, zvn, zwn     , tsn, tsa, jpts, nn_cen_h, nn_cen_v ) 
     135      CASE ( np_FCT )                                    ! FCT scheme      : 2nd / 4th order 
     136         CALL tra_adv_fct    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts, nn_fct_h, nn_fct_v ) 
     137      CASE ( np_FCT_zts )                                ! 2nd order FCT with vertical time-splitting 
     138         CALL tra_adv_fct_zts( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts        , nn_fct_zts ) 
     139      CASE ( np_MUS )                                    ! MUSCL 
     140         CALL tra_adv_mus    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb,      tsa, jpts        , ln_mus_ups )  
     141      CASE ( np_UBS )                                    ! UBS 
     142         CALL tra_adv_ubs    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts        , nn_ubs_v   ) 
     143      CASE ( np_QCK )                                    ! QUICKEST 
     144         CALL tra_adv_qck    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts                     ) 
     145      ! 
    151146      END SELECT 
    152147      ! 
    … …  
    157152      IF( nn_timing == 1 )  CALL timing_stop( 'tra_adv' ) 
    158153      ! 
    159       CALL wrk_dealloc( jpi, jpj, jpk, zun, zvn, zwn ) 
     154      CALL wrk_dealloc( jpi,jpj,jpk,  zun, zvn, zwn ) 
    160155      !                                           
    161156   END SUBROUTINE tra_adv 
    … …  
    169164      !!              tracer advection schemes and set nadv 
    170165      !!---------------------------------------------------------------------- 
    171       INTEGER ::   ioptio 
    172       INTEGER ::   ios                 ! Local integer output status for namelist read 
    173       !! 
    174       NAMELIST/namtra_adv/ ln_traadv_cen2 , ln_traadv_tvd,     & 
    175          &                 ln_traadv_muscl, ln_traadv_muscl2,  & 
    176          &                 ln_traadv_ubs  , ln_traadv_qck,     & 
    177          &                 ln_traadv_msc_ups, ln_traadv_tvd_zts 
    178       !!---------------------------------------------------------------------- 
    179  
    180       REWIND( numnam_ref )              ! Namelist namtra_adv in reference namelist : Tracer advection scheme 
     166      INTEGER ::   ioptio, ios   ! Local integers 
     167      ! 
     168      NAMELIST/namtra_adv/ ln_traadv_cen, nn_cen_h, nn_cen_v,               &   ! CEN 
     169         &                 ln_traadv_fct, nn_fct_h, nn_fct_v, nn_fct_zts,   &   ! FCT 
     170         &                 ln_traadv_mus,                     ln_mus_ups,   &   ! MUSCL 
     171         &                 ln_traadv_ubs,           nn_ubs_v,               &   ! UBS 
     172         &                 ln_traadv_qck                                        ! QCK 
     173      !!---------------------------------------------------------------------- 
     174      ! 
     175      !                                !==  Namelist  ==! 
     176      ! 
     177      REWIND( numnam_ref )                   ! Namelist namtra_adv in reference namelist : Tracer advection scheme 
    181178      READ  ( numnam_ref, namtra_adv, IOSTAT = ios, ERR = 901) 
    182179901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_adv in reference namelist', lwp ) 
    183  
    184       REWIND( numnam_cfg )              ! Namelist namtra_adv in configuration namelist : Tracer advection scheme 
     180      ! 
     181      REWIND( numnam_cfg )                   ! Namelist namtra_adv in configuration namelist : Tracer advection scheme 
    185182      READ  ( numnam_cfg, namtra_adv, IOSTAT = ios, ERR = 902 ) 
    186183902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_adv in configuration namelist', lwp ) 
    187184      IF(lwm) WRITE ( numond, namtra_adv ) 
    188185 
    189       IF(lwp) THEN                    ! Namelist print 
     186      IF(lwp) THEN                           ! Namelist print 
    190187         WRITE(numout,*) 
    191188         WRITE(numout,*) 'tra_adv_init : choice/control of the tracer advection scheme' 
    192189         WRITE(numout,*) '~~~~~~~~~~~' 
    193190         WRITE(numout,*) '   Namelist namtra_adv : chose a advection scheme for tracers' 
    194          WRITE(numout,*) '      2nd order advection scheme     ln_traadv_cen2    = ', ln_traadv_cen2 
    195          WRITE(numout,*) '      TVD advection scheme           ln_traadv_tvd     = ', ln_traadv_tvd 
    196          WRITE(numout,*) '      MUSCL  advection scheme        ln_traadv_muscl   = ', ln_traadv_muscl 
    197          WRITE(numout,*) '      MUSCL2 advection scheme        ln_traadv_muscl2  = ', ln_traadv_muscl2 
    198          WRITE(numout,*) '      UBS    advection scheme        ln_traadv_ubs     = ', ln_traadv_ubs 
    199          WRITE(numout,*) '      QUICKEST advection scheme      ln_traadv_qck     = ', ln_traadv_qck 
    200          WRITE(numout,*) '      upstream scheme within muscl   ln_traadv_msc_ups = ', ln_traadv_msc_ups 
    201          WRITE(numout,*) '      TVD advection scheme with zts  ln_traadv_tvd_zts = ', ln_traadv_tvd_zts 
    202       ENDIF 
    203  
    204       ioptio = 0                      ! Parameter control 
    205       IF( ln_traadv_cen2   )   ioptio = ioptio + 1 
    206       IF( ln_traadv_tvd    )   ioptio = ioptio + 1 
    207       IF( ln_traadv_muscl  )   ioptio = ioptio + 1 
    208       IF( ln_traadv_muscl2 )   ioptio = ioptio + 1 
    209       IF( ln_traadv_ubs    )   ioptio = ioptio + 1 
    210       IF( ln_traadv_qck    )   ioptio = ioptio + 1 
    211       IF( ln_traadv_tvd_zts)   ioptio = ioptio + 1 
    212       IF( lk_esopa         )   ioptio =          1 
    213  
    214       IF( ( ln_traadv_muscl .OR. ln_traadv_muscl2 .OR. ln_traadv_ubs .OR. ln_traadv_qck .OR. ln_traadv_tvd_zts )   & 
    215          .AND. ln_isfcav )   CALL ctl_stop( 'Only traadv_cen2 and traadv_tvd is compatible with ice shelf cavity') 
    216  
    217       IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE advection scheme in namelist namtra_adv' ) 
    218  
    219       !                              ! Set nadv 
    220       IF( ln_traadv_cen2   )   nadv =  1 
    221       IF( ln_traadv_tvd    )   nadv =  2 
    222       IF( ln_traadv_muscl  )   nadv =  3 
    223       IF( ln_traadv_muscl2 )   nadv =  4 
    224       IF( ln_traadv_ubs    )   nadv =  5 
    225       IF( ln_traadv_qck    )   nadv =  6 
    226       IF( ln_traadv_tvd_zts)   nadv =  7 
    227       IF( lk_esopa         )   nadv = -1 
    228  
    229       IF(lwp) THEN                   ! Print the choice 
     191         WRITE(numout,*) '      centered scheme                           ln_traadv_cen = ', ln_traadv_cen 
     192         WRITE(numout,*) '            horizontal 2nd/4th order               nn_cen_h   = ', nn_fct_h 
     193         WRITE(numout,*) '            vertical   2nd/4th order               nn_cen_v   = ', nn_fct_v 
     194         WRITE(numout,*) '      Flux Corrected Transport scheme           ln_traadv_fct = ', ln_traadv_fct 
     195         WRITE(numout,*) '            horizontal 2nd/4th order               nn_fct_h   = ', nn_fct_h 
     196         WRITE(numout,*) '            vertical   2nd/4th order               nn_fct_v   = ', nn_fct_v 
     197         WRITE(numout,*) '            2nd order + vertical sub-timestepping  nn_fct_zts = ', nn_fct_zts 
     198         WRITE(numout,*) '      MUSCL scheme                              ln_traadv_mus = ', ln_traadv_mus 
     199         WRITE(numout,*) '            + upstream scheme near river mouths    ln_mus_ups = ', ln_mus_ups 
     200         WRITE(numout,*) '      UBS scheme                                ln_traadv_ubs = ', ln_traadv_ubs 
     201         WRITE(numout,*) '            vertical   2nd/4th order               nn_ubs_v   = ', nn_ubs_v 
     202         WRITE(numout,*) '      QUICKEST scheme                           ln_traadv_qck = ', ln_traadv_qck 
     203      ENDIF 
     204 
     205      ioptio = 0                       !==  Parameter control  ==! 
     206      IF( ln_traadv_cen )   ioptio = ioptio + 1 
     207      IF( ln_traadv_fct )   ioptio = ioptio + 1 
     208      IF( ln_traadv_mus )   ioptio = ioptio + 1 
     209      IF( ln_traadv_ubs )   ioptio = ioptio + 1 
     210      IF( ln_traadv_qck )   ioptio = ioptio + 1 
     211      ! 
     212      IF( ioptio == 0 ) THEN 
     213         nadv = np_NO_adv 
     214         CALL ctl_warn( 'tra_adv_init: You are running without tracer advection.' ) 
     215      ENDIF 
     216      IF( ioptio /= 1 )   CALL ctl_stop( 'tra_adv_init: Choose ONE advection scheme in namelist namtra_adv' ) 
     217      ! 
     218      IF( ln_traadv_cen .AND. ( nn_cen_h /= 2 .AND. nn_cen_h /= 4 )   &          ! Centered 
     219                        .AND. ( nn_cen_v /= 2 .AND. nn_cen_v /= 4 )   ) THEN 
     220        CALL ctl_stop( 'tra_adv_init: CEN scheme, choose 2nd or 4th order' ) 
     221      ENDIF 
     222      IF( ln_traadv_fct .AND. ( nn_fct_h /= 2 .AND. nn_fct_h /= 4 )   &          ! FCT 
     223                        .AND. ( nn_fct_v /= 2 .AND. nn_fct_v /= 4 )   ) THEN 
     224        CALL ctl_stop( 'tra_adv_init: FCT scheme, choose 2nd or 4th order' ) 
     225      ENDIF 
     226      IF( ln_traadv_fct .AND. nn_fct_zts > 0 ) THEN 
     227         IF( nn_fct_h == 4 ) THEN 
     228            nn_fct_h = 2 
     229            CALL ctl_stop( 'tra_adv_init: force 2nd order FCT scheme, 4th order does not exist with sub-timestepping' ) 
     230         ENDIF 
     231         IF( lk_vvl ) THEN 
     232            CALL ctl_stop( 'tra_adv_init: vertical sub-timestepping not allow in non-linear free surface' ) 
     233         ENDIF 
     234         IF( nn_fct_zts == 1 )   CALL ctl_warn( 'tra_adv_init: FCT with ONE sub-timestep = FCT without sub-timestep' ) 
     235      ENDIF 
     236      IF( ln_traadv_ubs .AND. ( nn_ubs_v /= 2 .AND. nn_ubs_v /= 4 )   ) THEN     ! UBS 
     237        CALL ctl_stop( 'tra_adv_init: UBS scheme, choose 2nd or 4th order' ) 
     238      ENDIF 
     239      IF( ln_traadv_ubs .AND. nn_ubs_v == 4 ) THEN 
     240         CALL ctl_warn( 'tra_adv_init: UBS scheme, only 2nd FCT scheme available on the vertical. It will be used' ) 
     241      ENDIF 
     242      IF( ln_isfcav ) THEN                                                       ! ice-shelf cavities 
     243         IF(  ln_traadv_cen .AND. nn_cen_v /= 4    .OR.   &                            ! NO 4th order with ISF 
     244            & ln_traadv_fct .AND. nn_fct_v /= 4   )   CALL ctl_stop( 'tra_adv_init: 4th order COMPACT scheme not allowed with ISF' ) 
     245      ENDIF 
     246      ! 
     247      !                                !==  used advection scheme  ==!   
     248      !                                      ! set nadv 
     249      IF( ln_traadv_cen                      )   nadv = np_CEN 
     250      IF( ln_traadv_fct                      )   nadv = np_FCT 
     251      IF( ln_traadv_fct .AND. nn_fct_zts > 0 )   nadv = np_FCT_zts 
     252      IF( ln_traadv_mus                      )   nadv = np_MUS 
     253      IF( ln_traadv_ubs                      )   nadv = np_UBS 
     254      IF( ln_traadv_qck                      )   nadv = np_QCK 
     255 
     256      IF(lwp) THEN                           ! Print the choice 
    230257         WRITE(numout,*) 
    231          IF( nadv ==  1 )   WRITE(numout,*) '         2nd order scheme is used' 
    232          IF( nadv ==  2 )   WRITE(numout,*) '         TVD       scheme is used' 
    233          IF( nadv ==  3 )   WRITE(numout,*) '         MUSCL     scheme is used' 
    234          IF( nadv ==  4 )   WRITE(numout,*) '         MUSCL2    scheme is used' 
    235          IF( nadv ==  5 )   WRITE(numout,*) '         UBS       scheme is used' 
    236          IF( nadv ==  6 )   WRITE(numout,*) '         QUICKEST  scheme is used' 
    237          IF( nadv ==  7 )   WRITE(numout,*) '         TVD ZTS   scheme is used' 
    238          IF( nadv == -1 )   WRITE(numout,*) '         esopa test: use all advection scheme' 
    239       ENDIF 
    240       ! 
    241       CALL tra_adv_mle_init          ! initialisation of the Mixed Layer Eddy parametrisation (MLE) 
     258         IF( nadv == np_NO_adv  )   WRITE(numout,*) '         NO T-S advection' 
     259         IF( nadv == np_CEN     )   WRITE(numout,*) '         CEN      scheme is used. Horizontal order: ', nn_cen_h,   & 
     260            &                                                                        ' Vertical   order: ', nn_cen_v 
     261         IF( nadv == np_FCT     )   WRITE(numout,*) '         FCT      scheme is used. Horizontal order: ', nn_fct_h,   & 
     262            &                                                                        ' Vertical   order: ', nn_fct_v 
     263         IF( nadv == np_FCT_zts )   WRITE(numout,*) '         use 2nd order FCT with ', nn_fct_zts,'vertical sub-timestepping' 
     264         IF( nadv == np_MUS     )   WRITE(numout,*) '         MUSCL    scheme is used' 
     265         IF( nadv == np_UBS     )   WRITE(numout,*) '         UBS      scheme is used' 
     266         IF( nadv == np_QCK     )   WRITE(numout,*) '         QUICKEST scheme is used' 
     267      ENDIF 
     268      ! 
     269      CALL tra_adv_mle_init            !== initialisation of the Mixed Layer Eddy parametrisation (MLE)  ==! 
    242270      ! 
    243271   END SUBROUTINE tra_adv_init 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mle.F90

    r5215 r5836  
    2828   PUBLIC   tra_adv_mle_init   ! routine called in traadv.F90 
    2929 
    30    !                                               !!* namelist namtra_adv_mle * 
     30   !                                       !!* namelist namtra_adv_mle * 
    3131   LOGICAL, PUBLIC ::   ln_mle              ! flag to activate the Mixed Layer Eddy (MLE) parameterisation 
    3232   INTEGER         ::   nn_mle              ! MLE type: =0 standard Fox-Kemper ; =1 new formulation 
    … …  
    3434   INTEGER         ::   nn_conv             ! =1 no MLE in case of convection ; =0 always MLE 
    3535   REAL(wp)        ::   rn_ce               ! MLE coefficient 
    36    !                                           ! parameters used in nn_mle = 0 case 
     36   !                                        ! parameters used in nn_mle = 0 case 
    3737   REAL(wp)        ::   rn_lf                  ! typical scale of mixed layer front 
    38    REAL(wp)        ::   rn_time             ! time scale for mixing momentum across the mixed layer 
    39    !                                             ! parameters used in nn_mle = 1 case 
    40    REAL(wp)        ::   rn_lat                   ! reference latitude for a 5 km scale of ML front 
    41    REAL(wp)        ::   rn_rho_c_mle         ! Density criterion for definition of MLD used by FK 
     38   REAL(wp)        ::   rn_time                ! time scale for mixing momentum across the mixed layer 
     39   !                                        ! parameters used in nn_mle = 1 case 
     40   REAL(wp)        ::   rn_lat                 ! reference latitude for a 5 km scale of ML front 
     41   REAL(wp)        ::   rn_rho_c_mle           ! Density criterion for definition of MLD used by FK 
    4242 
    4343   REAL(wp) ::   r5_21 = 5.e0 / 21.e0   ! factor used in mle streamfunction computation 
    … …  
    5252#  include "vectopt_loop_substitute.h90" 
    5353   !!---------------------------------------------------------------------- 
    54    !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
     54   !! NEMO/OPA 4.0 , NEMO Consortium (2015) 
    5555   !! $Id$ 
    5656   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    … …  
    8080      !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 
    8181      !!---------------------------------------------------------------------- 
    82       ! 
    8382      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index 
    8483      INTEGER                         , INTENT(in   ) ::   kit000     ! first time step index 
    … …  
    9392      REAL(wp) ::   zcvw, zmvw   !   -      - 
    9493      REAL(wp) ::   zc                                     !   -      - 
    95  
     94      ! 
    9695      INTEGER  ::   ii, ij, ik              ! local integers 
    9796      INTEGER, DIMENSION(3) ::   ilocu      ! 
    … …  
    101100      INTEGER, POINTER, DIMENSION(:,:) :: inml_mle 
    102101      !!---------------------------------------------------------------------- 
    103  
     102      ! 
    104103      IF( nn_timing == 1 )  CALL timing_start('tra_adv_mle') 
    105104      CALL wrk_alloc( jpi, jpj, zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH) 
    … …  
    171170         DO jj = 1, jpjm1 
    172171            DO ji = 1, fs_jpim1   ! vector opt. 
    173                zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2u(ji,jj)                                            & 
    174                   &           * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp          , e1u(ji,jj)                )   & 
    175                   &           / (         e1u(ji,jj)          * MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) )   ) 
     172               zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2_e1u(ji,jj)                                            & 
     173                  &           * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )   & 
     174                  &           / (  MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) )   ) 
    176175                  ! 
    177                zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj)  * e1v(ji,jj)                                            & 
    178                   &           * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp          , e2v(ji,jj)                )   & 
    179                   &           / (         e2v(ji,jj)          * MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) )   ) 
     176               zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj)  * e1_e2v(ji,jj)                                            & 
     177                  &           * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )   & 
     178                  &           / (  MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) )   ) 
    180179            END DO 
    181180         END DO 
    … …  
    184183         DO jj = 1, jpjm1 
    185184            DO ji = 1, fs_jpim1   ! vector opt. 
    186                zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2u(ji,jj) / e1u(ji,jj)          & 
     185               zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2_e1u(ji,jj)               & 
    187186                  &                  * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) 
    188187                  ! 
    189                zpsim_v(ji,jj) = rc_f *   zhv(ji,jj)   * zhv(ji,jj)   * e1v(ji,jj) / e2v(ji,jj)          & 
     188               zpsim_v(ji,jj) = rc_f *   zhv(ji,jj)   * zhv(ji,jj)   * e1_e2v(ji,jj)               & 
    190189                  &                  * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) 
    191190            END DO 
    … …  
    252251         ! divide by cross distance to give streamfunction with dimensions m^2/s 
    253252         DO jk = 1, ikmax+1 
    254             zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk)/e2u(:,:) 
    255             zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk)/e1v(:,:) 
     253            zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk) * r1_e2u(:,:) 
     254            zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk) * r1_e1v(:,:) 
    256255         END DO 
    257256         CALL iom_put( "psiu_mle", zpsi_uw )    ! i-mle streamfunction 
    … …  
    281280      NAMELIST/namtra_adv_mle/ ln_mle , nn_mle, rn_ce, rn_lf, rn_time, rn_lat, nn_mld_uv, nn_conv, rn_rho_c_mle 
    282281      !!---------------------------------------------------------------------- 
    283  
    284282 
    285283      REWIND( numnam_ref )              ! Namelist namtra_adv_mle in reference namelist : Tracer advection scheme 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90

    r5147 r5836  
    102102         IF(lwp) WRITE(numout,*) 
    103103      ENDIF 
     104      ! 
    104105      l_trd = .FALSE. 
    105106      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE. 
    … …  
    130131      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend  
    131132      !! 
    132       INTEGER  :: ji, jj, jk, jn   ! dummy loop indices 
    133       REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk   ! local scalars 
    134       REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zfu, zfc, zfd 
     133      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
     134      REAL(wp) ::   ztra, zbtr, zdir, zdx, zdt, zmsk   ! local scalars 
     135      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwx, zfu, zfc, zfd 
    135136      !---------------------------------------------------------------------- 
    136137      ! 
    … …  
    139140      DO jn = 1, kjpt                                            ! tracer loop 
    140141         !                                                       ! =========== 
    141          zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0   
    142          zfd(:,:,:) = 0.0     ;   zwx(:,:,:) = 0.0      
    143          !                                                   
    144          DO jk = 1, jpkm1                                 
    145             !                                              
    146             !--- Computation of the ustream and downstream value of the tracer and the mask 
     142         zfu(:,:,:) = 0._wp     ;   zfc(:,:,:) = 0._wp  
     143         zfd(:,:,:) = 0._wp     ;   zwx(:,:,:) = 0._wp    
     144         ! 
     145!!gm why not using a SHIFT instruction... 
     146         DO jk = 1, jpkm1     !--- Computation of the ustream and downstream value of the tracer and the mask 
    147147            DO jj = 2, jpjm1 
    148148               DO ji = fs_2, fs_jpim1   ! vector opt. 
    149                   ! Upstream in the x-direction for the tracer 
    150                   zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn) 
    151                   ! Downstream in the x-direction for the tracer 
    152                   zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn) 
     149                  zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn)        ! Upstream   in the x-direction for the tracer 
     150                  zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn)        ! Downstream in the x-direction for the tracer 
    153151               END DO 
    154152            END DO 
    … …  
    159157         ! Horizontal advective fluxes 
    160158         ! --------------------------- 
    161          ! 
    162159         DO jk = 1, jpkm1                              
    163160            DO jj = 2, jpjm1 
    … …  
    220217            DO jj = 2, jpjm1 
    221218               DO ji = fs_2, fs_jpim1   ! vector opt.   
    222                   zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     219                  zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    223220                  ! horizontal advective trends 
    224221                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) ) 
    … …  
    344341            DO jj = 2, jpjm1 
    345342               DO ji = fs_2, fs_jpim1   ! vector opt.   
    346                   zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     343                  zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    347344                  ! horizontal advective trends 
    348345                  ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) ) 
    … …  
    380377      ! 
    381378      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    382       REAL(wp) ::   zbtr , ztra      ! local scalars 
    383379      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz 
    384380      !!---------------------------------------------------------------------- 
    385381      ! 
    386       CALL wrk_alloc( jpi, jpj, jpk, zwz ) 
     382      CALL wrk_alloc( jpi,jpj,jpk,   zwz ) 
     383      ! 
     384      !                          ! surface & bottom values  
     385      IF( lk_vvl )   zwz(:,:, 1 ) = 0._wp             ! set to zero one for all 
     386                     zwz(:,:,jpk) = 0._wp             ! except at the surface in linear free surface 
     387      ! 
    387388      !                                                          ! =========== 
    388389      DO jn = 1, kjpt                                            ! tracer loop 
    389390         !                                                       ! =========== 
    390          ! 1. Bottom value : flux set to zero 
    391          zwz(:,:,jpk) = 0.e0             ! Bottom value : flux set to zero 
    392          ! 
    393          !                                 ! Surface value 
    394          IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0.e0                      ! Variable volume : flux set to zero 
    395          ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn)   ! Constant volume : advective flux through the surface 
     391         ! 
     392         DO jk = 2, jpkm1                    !* Interior point   (w-masked 2nd order centered flux) 
     393            DO jj = 2, jpjm1 
     394               DO ji = fs_2, fs_jpim1   ! vector opt. 
     395                  zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) ) * wmask(ji,jj,jk) 
     396               END DO 
     397            END DO 
     398         END DO 
     399         IF(.NOT.lk_vvl ) THEN               !* top value   (only in linear free surf. as zwz is multiplied by wmask) 
     400            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean) 
     401               DO jj = 1, jpj 
     402                  DO ji = 1, jpi 
     403                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
     404                  END DO 
     405               END DO    
     406            ELSE                                   ! no ice-shelf cavities (only ocean surface) 
     407               zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn) 
     408            ENDIF 
    396409         ENDIF 
    397410         ! 
    398          DO jk = 2, jpkm1                  ! Interior point: second order centered tracer flux at w-point 
     411         DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==! 
    399412            DO jj = 2, jpjm1 
    400413               DO ji = fs_2, fs_jpim1   ! vector opt. 
    401                   zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) ) 
    402                END DO 
    403             END DO 
    404          END DO 
    405          ! 
    406          DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==! 
    407             DO jj = 2, jpjm1 
    408                DO ji = fs_2, fs_jpim1   ! vector opt. 
    409                   zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    410                   ! k- vertical advective trends  
    411                   ztra = - zbtr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )  
    412                   ! added to the general tracer trends 
    413                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 
     414                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   & 
     415                     &                                / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    414416               END DO 
    415417            END DO 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90

    r5147 r5836  
    1616   USE trc_oce        ! share passive tracers/Ocean variables 
    1717   USE trd_oce        ! trends: ocean variables 
     18   USE traadv_fct      ! acces to routine interp_4th_cpt  
    1819   USE trdtra         ! trends manager: tracers  
    1920   USE dynspg_oce     ! choice/control of key cpp for surface pressure gradient 
    … …  
    3839#  include "vectopt_loop_substitute.h90" 
    3940   !!---------------------------------------------------------------------- 
    40    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     41   !! NEMO/OPA 3.7 , NEMO Consortium (2015) 
    4142   !! $Id$ 
    4243   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    4445CONTAINS 
    4546 
    46    SUBROUTINE tra_adv_ubs ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      & 
    47       &                                       ptb, ptn, pta, kjpt ) 
     47   SUBROUTINE tra_adv_ubs( kt, kit000, cdtype, p2dt, pun, pvn, pwn,          & 
     48      &                                                ptb, ptn, pta, kjpt, kn_ubs_v ) 
    4849      !!---------------------------------------------------------------------- 
    4950      !!                  ***  ROUTINE tra_adv_ubs  *** 
    … …  
    5253      !!      and add it to the general trend of passive tracer equations. 
    5354      !! 
    54       !! ** Method  :   The upstream biased scheme (UBS) is based on a 3rd order 
     55      !! ** Method  :   The 3rd order Upstream Biased Scheme (UBS) is based on an 
    5556      !!      upstream-biased parabolic interpolation (Shchepetkin and McWilliams 2005) 
    5657      !!      It is only used in the horizontal direction. 
    … …  
    6162      !!      where zltu is the second derivative of the before temperature field: 
    6263      !!          zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ] 
    63       !!      This results in a dissipatively dominant (i.e. hyper-diffusive)  
     64      !!        This results in a dissipatively dominant (i.e. hyper-diffusive)  
    6465      !!      truncation error. The overall performance of the advection scheme  
    6566      !!      is similar to that reported in (Farrow and Stevens, 1995).  
    66       !!      For stability reasons, the first term of the fluxes which corresponds 
     67      !!        For stability reasons, the first term of the fluxes which corresponds 
    6768      !!      to a second order centered scheme is evaluated using the now velocity  
    6869      !!      (centered in time) while the second term which is the diffusive part  
    6970      !!      of the scheme, is evaluated using the before velocity (forward in time).  
    7071      !!      Note that UBS is not positive. Do not use it on passive tracers. 
    71       !!                On the vertical, the advection is evaluated using a TVD scheme, 
    72       !!      as the UBS have been found to be too diffusive. 
     72      !!                On the vertical, the advection is evaluated using a FCT scheme, 
     73      !!      as the UBS have been found to be too diffusive.  
     74!!gm  !!                kn_ubs_v argument (not coded for the moment) 
     75      !!      controles whether the FCT is based on a 2nd order centrered scheme (kn_ubs_v=2)  
     76      !!      or on a 4th order compact scheme (kn_ubs_v=4). 
    7377      !! 
    7478      !! ** Action : - update (pta) with the now advective tracer trends 
    … …  
    8185      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator) 
    8286      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers 
     87      INTEGER                              , INTENT(in   ) ::   kn_ubs_v        ! number of tracers 
    8388      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step 
    8489      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean transport components 
    … …  
    95100      IF( nn_timing == 1 )  CALL timing_start('tra_adv_ubs') 
    96101      ! 
    97       CALL wrk_alloc( jpi, jpj, jpk, ztu, ztv, zltu, zltv, zti, ztw ) 
     102      CALL wrk_alloc( jpi,jpj,jpk,  ztu, ztv, zltu, zltv, zti, ztw ) 
    98103      ! 
    99104      IF( kt == kit000 )  THEN 
    … …  
    106111      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 
    107112      ! 
     113      zltu(:,:,jpk) = 0._wp   ;   zltv(:,:,jpk) = 0._wp     ! Bottom value : set to zero one for all 
     114      ztw (:,:,jpk) = 0._wp   ;   zti (:,:,jpk) = 0._wp 
     115      IF( lk_vvl )   ztw(:,:, 1 ) = 0._wp                   ! surface value: set to zero only in vvl case 
     116      ! 
    108117      !                                                          ! =========== 
    109118      DO jn = 1, kjpt                                            ! tracer loop 
    110119         !                                                       ! =========== 
    111          ! 1. Bottom value : flux set to zero 
    112          ! ---------------------------------- 
    113          zltu(:,:,jpk) = 0.e0       ;      zltv(:,:,jpk) = 0.e0 
    114120         !                                               
    115          DO jk = 1, jpkm1                                 ! Horizontal slab 
    116             !                                    
    117             !  Laplacian 
    118             DO jj = 1, jpjm1            ! First derivative (gradient) 
     121         DO jk = 1, jpkm1        !==  horizontal laplacian of before tracer ==! 
     122            DO jj = 1, jpjm1              ! First derivative (masked gradient) 
    119123               DO ji = 1, fs_jpim1   ! vector opt. 
    120                   zeeu = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk) 
    121                   zeev = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk) 
     124                  zeeu = e2_e1u(ji,jj) * fse3u(ji,jj,jk) * umask(ji,jj,jk) 
     125                  zeev = e1_e2v(ji,jj) * fse3v(ji,jj,jk) * vmask(ji,jj,jk) 
    122126                  ztu(ji,jj,jk) = zeeu * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) 
    123127                  ztv(ji,jj,jk) = zeev * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) 
    124128               END DO 
    125129            END DO 
    126             DO jj = 2, jpjm1            ! Second derivative (divergence) 
     130            DO jj = 2, jpjm1              ! Second derivative (divergence) 
    127131               DO ji = fs_2, fs_jpim1   ! vector opt. 
    128                   zcoef = 1. / ( 6. * fse3t(ji,jj,jk) ) 
     132                  zcoef = 1._wp / ( 6._wp * fse3t(ji,jj,jk) ) 
    129133                  zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef 
    130134                  zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef 
    … …  
    132136            END DO 
    133137            !                                     
    134          END DO                                           ! End of slab          
     138         END DO          
    135139         CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn) 
    136  
    137140         !     
    138          !  Horizontal advective fluxes                
    139          DO jk = 1, jpkm1                                 ! Horizontal slab 
     141         DO jk = 1, jpkm1        !==  Horizontal advective fluxes  ==!     (UBS) 
    140142            DO jj = 1, jpjm1 
    141143               DO ji = 1, fs_jpim1   ! vector opt. 
    142                   ! upstream transport (x2) 
    143                   zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) 
     144                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )      ! upstream transport (x2) 
    144145                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) 
    145146                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) 
    146147                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) 
    147                   ! 2nd order centered advective fluxes (x2) 
     148                  !                                                  ! 2nd order centered advective fluxes (x2) 
    148149                  zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn) ) 
    149150                  zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) ) 
    150                   ! UBS advective fluxes 
     151                  !                                                  ! UBS advective fluxes 
    151152                  ztu(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) ) 
    152153                  ztv(ji,jj,jk) = 0.5 * ( zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) ) 
    153154               END DO 
    154155            END DO 
    155          END DO                                           ! End of slab          
    156  
    157          zltu(:,:,:) = pta(:,:,:,jn)      ! store pta trends 
    158  
    159          DO jk = 1, jpkm1                 ! Horizontal advective trends 
     156         END DO          
     157         ! 
     158         zltu(:,:,:) = pta(:,:,:,jn)      ! store the initial trends before its update 
     159         ! 
     160         DO jk = 1, jpkm1        !==  add the horizontal advective trend  ==! 
    160161            DO jj = 2, jpjm1 
    161162               DO ji = fs_2, fs_jpim1   ! vector opt. 
    … …  
    166167            END DO 
    167168            !                                              
    168          END DO                                           !   End of slab 
    169  
    170          ! Horizontal trend used in tra_adv_ztvd subroutine 
    171          zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:) 
    172  
     169         END DO 
     170         ! 
     171         zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:)    ! Horizontal advective trend used in vertical 2nd order FCT case 
     172         !                                            ! and/or in trend diagnostic (l_trd=T)  
    173173         !                 
    174174         IF( l_trd ) THEN                  ! trend diagnostics 
    … …  
    181181            IF( jn == jp_sal )  str_adv(:) = ptr_sj( ztv(:,:,:) ) 
    182182         ENDIF 
    183           
    184          ! TVD scheme for the vertical direction   
    185          ! ---------------------- 
    186          IF( l_trd )   zltv(:,:,:) = pta(:,:,:,jn)          ! store pta if trend diag. 
    187  
    188          !  Bottom value : flux set to zero 
    189          ztw(:,:,jpk) = 0.e0   ;   zti(:,:,jpk) = 0.e0 
    190  
    191          ! Surface value 
    192          IF( lk_vvl ) THEN   ;   ztw(:,:,1) = 0.e0                      ! variable volume : flux set to zero 
    193          ELSE                ;   ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)   ! free constant surface  
    194          ENDIF 
    195          !  upstream advection with initial mass fluxes & intermediate update 
    196          ! ------------------------------------------------------------------- 
    197          ! Interior value 
    198          DO jk = 2, jpkm1 
    199             DO jj = 1, jpj 
    200                DO ji = 1, jpi 
    201                    zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 
    202                    zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 
    203                    ztw(ji,jj,jk) = 0.5 * (  zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn)  ) 
    204                END DO 
    205             END DO 
    206          END DO  
    207          ! update and guess with monotonic sheme 
    208          DO jk = 1, jpkm1 
    209             z2dtt = p2dt(jk) 
    210             DO jj = 2, jpjm1 
    211                DO ji = fs_2, fs_jpim1   ! vector opt. 
    212                   zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    213                   ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr 
    214                   pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn) +  ztak  
    215                   zti(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk) 
    216                END DO 
    217             END DO 
    218          END DO 
    219          ! 
    220          CALL lbc_lnk( zti, 'T', 1. )      ! Lateral boundary conditions on zti, zsi   (unchanged sign) 
    221  
    222          !  antidiffusive flux : high order minus low order 
    223          ztw(:,:,1) = 0.e0       ! Surface value 
    224          DO jk = 2, jpkm1        ! Interior value 
    225             DO jj = 1, jpj 
    226                DO ji = 1, jpi 
    227                   ztw(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - ztw(ji,jj,jk) 
    228                END DO 
    229             END DO 
    230          END DO 
    231          ! 
    232          CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt )      !  monotonicity algorithm 
    233  
    234          !  final trend with corrected fluxes 
    235          DO jk = 1, jpkm1 
     183         ! 
     184         !                       !== vertical advective trend  ==! 
     185         ! 
     186         SELECT CASE( kn_ubs_v )       ! select the vertical advection scheme 
     187         ! 
     188         CASE(  2  )                   ! 2nd order FCT  
     189            !          
     190            IF( l_trd )   zltv(:,:,:) = pta(:,:,:,jn)          ! store pta if trend diag. 
     191            ! 
     192            !                          !*  upstream advection with initial mass fluxes & intermediate update  ==! 
     193            DO jk = 2, jpkm1                 ! Interior value (w-masked) 
     194               DO jj = 1, jpj 
     195                  DO ji = 1, jpi 
     196                     zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 
     197                     zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 
     198                     ztw(ji,jj,jk) = 0.5_wp * (  zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn)  ) * wmask(ji,jj,jk) 
     199                  END DO 
     200               END DO 
     201            END DO  
     202            IF(.NOT.lk_vvl ) THEN            ! top ocean value (only in linear free surface as ztw has been w-masked) 
     203               IF( ln_isfcav ) THEN                ! top of the ice-shelf cavities and at the ocean surface 
     204                  DO jj = 1, jpj 
     205                     DO ji = 1, jpi 
     206                        ztw(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
     207                     END DO 
     208                  END DO    
     209               ELSE                                ! no cavities: only at the ocean surface 
     210                  ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 
     211               ENDIF 
     212            ENDIF 
     213            ! 
     214            DO jk = 1, jpkm1           !* trend and after field with monotonic scheme 
     215               z2dtt = p2dt(jk) 
     216               DO jj = 2, jpjm1 
     217                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     218                     ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 
     219                     pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn) +  ztak  
     220                     zti(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk) 
     221                  END DO 
     222               END DO 
     223            END DO 
     224            CALL lbc_lnk( zti, 'T', 1. )      ! Lateral boundary conditions on zti, zsi   (unchanged sign) 
     225            ! 
     226            !                          !*  anti-diffusive flux : high order minus low order 
     227            DO jk = 2, jpkm1        ! Interior value  (w-masked) 
     228               DO jj = 1, jpj 
     229                  DO ji = 1, jpi 
     230                     ztw(ji,jj,jk) = (   0.5_wp * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )   & 
     231                        &              - ztw(ji,jj,jk)   ) * wmask(ji,jj,jk) 
     232                  END DO 
     233               END DO 
     234            END DO 
     235            !                                            ! top ocean value: high order == upstream  ==>>  zwz=0 
     236            IF(.NOT.lk_vvl )   ztw(:,:, 1 ) = 0._wp      ! only ocean surface as interior zwz values have been w-masked 
     237            ! 
     238            CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt )      !  monotonicity algorithm 
     239            ! 
     240         CASE(  4  )                               ! 4th order COMPACT 
     241            CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )         ! 4th order compact interpolation of T at w-point 
     242            DO jk = 2, jpkm1 
     243               DO jj = 2, jpjm1 
     244                  DO ji = fs_2, fs_jpim1 
     245                     ztw(ji,jj,jk) = pwn(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk) 
     246                  END DO 
     247               END DO 
     248            END DO 
     249            IF(.NOT.lk_vvl )   ztw(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn)     !!gm ISF & 4th COMPACT doesn't work 
     250            ! 
     251         END SELECT 
     252         ! 
     253         DO jk = 1, jpkm1        !  final trend with corrected fluxes 
    236254            DO jj = 2, jpjm1  
    237255               DO ji = fs_2, fs_jpim1   ! vector opt.    
    238                   zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    239                   ! k- vertical advective trends   
    240                   ztra = - zbtr * ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) 
    241                   ! added to the general tracer trends 
    242                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 
    243                END DO 
    244             END DO 
    245          END DO 
    246  
    247          !  Save the final vertical advective trends 
    248          IF( l_trd )  THEN                        ! vertical advective trend diagnostics 
     256                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 
     257               END DO 
     258            END DO 
     259         END DO 
     260         ! 
     261         IF( l_trd )  THEN       ! vertical advective trend diagnostics 
    249262            DO jk = 1, jpkm1                       ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w]) 
    250263               DO jj = 2, jpjm1 
    251264                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    252                      zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    253                      z_hdivn = (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  ) * zbtr 
    254                      zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk) + ptn(ji,jj,jk,jn) * z_hdivn 
     265                     zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk)                          & 
     266                        &           + ptn(ji,jj,jk,jn) * (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  )   & 
     267                        &                              / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 
    255268                  END DO 
    256269               END DO 
    … …  
    261274      END DO 
    262275      ! 
    263       CALL wrk_dealloc( jpi, jpj, jpk, ztu, ztv, zltu, zltv, zti, ztw ) 
     276      CALL wrk_dealloc( jpi,jpj,jpk,  ztu, ztv, zltu, zltv, zti, ztw ) 
    264277      ! 
    265278      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_ubs') 
    … …  
    294307      IF( nn_timing == 1 )  CALL timing_start('nonosc_z') 
    295308      ! 
    296       CALL wrk_alloc( jpi, jpj, jpk, zbetup, zbetdo ) 
     309      CALL wrk_alloc( jpi,jpj,jpk,  zbetup, zbetdo ) 
    297310      ! 
    298311      zbig  = 1.e+40_wp 
    299312      zrtrn = 1.e-15_wp 
    300313      zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp 
    301  
     314      ! 
    302315      ! Search local extrema 
    303316      ! -------------------- 
    304       ! large negative value (-zbig) inside land 
     317      !                    ! large negative value (-zbig) inside land 
    305318      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) ) 
    306319      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) ) 
    307       ! search maximum in neighbourhood 
    308       DO jk = 1, jpkm1 
     320      ! 
     321      DO jk = 1, jpkm1     ! search maximum in neighbourhood 
    309322         ikm1 = MAX(jk-1,1) 
    310323         DO jj = 2, jpjm1 
    … …  
    316329         END DO 
    317330      END DO 
    318       ! large positive value (+zbig) inside land 
     331      !                    ! large positive value (+zbig) inside land 
    319332      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) ) 
    320333      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) ) 
    321       ! search minimum in neighbourhood 
    322       DO jk = 1, jpkm1 
     334      ! 
     335      DO jk = 1, jpkm1     ! search minimum in neighbourhood 
    323336         ikm1 = MAX(jk-1,1) 
    324337         DO jj = 2, jpjm1 
    … …  
    330343         END DO 
    331344      END DO 
    332  
    333       ! restore masked values to zero 
     345      !                    ! restore masked values to zero 
    334346      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) 
    335347      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) 
    336  
    337  
    338       ! 2. Positive and negative part of fluxes and beta terms 
    339       ! ------------------------------------------------------ 
    340  
     348      ! 
     349      ! Positive and negative part of fluxes and beta terms 
     350      ! --------------------------------------------------- 
    341351      DO jk = 1, jpkm1 
    342352         z2dtt = p2dt(jk) 
    … …  
    347357               zneg = MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) ) 
    348358               ! up & down beta terms 
    349                zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt 
     359               zbt = e1e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt 
    350360               zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt 
    351361               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt 
    … …  
    353363         END DO 
    354364      END DO 
     365      ! 
    355366      ! monotonic flux in the k direction, i.e. pcc 
    356367      ! ------------------------------------------- 
    … …  
    366377      END DO 
    367378      ! 
    368       CALL wrk_dealloc( jpi, jpj, jpk, zbetup, zbetdo ) 
     379      CALL wrk_dealloc( jpi,jpj,jpk,  zbetup, zbetdo ) 
    369380      ! 
    370381      IF( nn_timing == 1 )  CALL timing_stop('nonosc_z') 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90

    r4990 r5836  
    1414   !!             -   ! 2013-04  (F. Roquet, G. Madec)  use of eosbn2 instead of local hard coded alpha and beta 
    1515   !!---------------------------------------------------------------------- 
    16 #if   defined key_trabbl   ||   defined key_esopa 
     16#if   defined key_trabbl 
    1717   !!---------------------------------------------------------------------- 
    1818   !!   'key_trabbl'   or                             bottom boundary layer 
    … …  
    2929   USE phycst         ! physical constant 
    3030   USE eosbn2         ! equation of state 
    31    USE trd_oce     ! trends: ocean variables 
     31   USE trd_oce        ! trends: ocean variables 
    3232   USE trdtra         ! trends: active tracers 
    3333   ! 
    … …  
    198198         DO jj = 1, jpj 
    199199            DO ji = 1, jpi 
    200                ik = mbkt(ji,jj)                              ! bottom T-level index 
    201                zptb(ji,jj) = ptb(ji,jj,ik,jn)       ! bottom before T and S 
     200               ik = mbkt(ji,jj)                             ! bottom T-level index 
     201               zptb(ji,jj) = ptb(ji,jj,ik,jn)               ! bottom before T and S 
    202202            END DO 
    203203         END DO 
    … …  
    205205         DO jj = 2, jpjm1                                    ! Compute the trend 
    206206            DO ji = 2, jpim1 
    207                ik = mbkt(ji,jj)                              ! bottom T-level index 
    208                zbtr = r1_e12t(ji,jj)  / fse3t(ji,jj,ik) 
    209                pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                         & 
    210                   &               + (   ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )   & 
    211                   &                   - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )   & 
    212                   &                   + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )   & 
    213                   &                   - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )   ) * zbtr 
     207               ik = mbkt(ji,jj)                            ! bottom T-level index 
     208               pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                  & 
     209                  &             + (  ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )     & 
     210                  &                - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )     & 
     211                  &                + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )     & 
     212                  &                - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )  )  & 
     213                  &             / ( e1e2t(ji,jj) * fse3t(ji,jj,ik) ) 
    214214            END DO 
    215215         END DO 
    … …  
    263263                  ! 
    264264                  !                                               ! up  -slope T-point (shelf bottom point) 
    265                   zbtr = r1_e12t(iis,jj) / fse3t(iis,jj,ikus) 
     265                  zbtr = r1_e1e2t(iis,jj) / fse3t(iis,jj,ikus) 
    266266                  ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr 
    267267                  pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra 
    268268                  ! 
    269269                  DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column) 
    270                      zbtr = r1_e12t(iid,jj) / fse3t(iid,jj,jk) 
     270                     zbtr = r1_e1e2t(iid,jj) / fse3t(iid,jj,jk) 
    271271                     ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr 
    272272                     pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra 
    273273                  END DO 
    274274                  ! 
    275                   zbtr = r1_e12t(iid,jj) / fse3t(iid,jj,ikud) 
     275                  zbtr = r1_e1e2t(iid,jj) / fse3t(iid,jj,ikud) 
    276276                  ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr 
    277277                  pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra 
    … …  
    285285                  ! 
    286286                  ! up  -slope T-point (shelf bottom point) 
    287                   zbtr = r1_e12t(ji,ijs) / fse3t(ji,ijs,ikvs) 
     287                  zbtr = r1_e1e2t(ji,ijs) / fse3t(ji,ijs,ikvs) 
    288288                  ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr 
    289289                  pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra 
    290290                  ! 
    291291                  DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column) 
    292                      zbtr = r1_e12t(ji,ijd) / fse3t(ji,ijd,jk) 
     292                     zbtr = r1_e1e2t(ji,ijd) / fse3t(ji,ijd,jk) 
    293293                     ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr 
    294294                     pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn)  + ztra 
    295295                  END DO 
    296296                  !                                               ! down-slope T-point (deep bottom point) 
    297                   zbtr = r1_e12t(ji,ijd) / fse3t(ji,ijd,ikvd) 
     297                  zbtr = r1_e1e2t(ji,ijd) / fse3t(ji,ijd,ikvd) 
    298298                  ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr 
    299299                  pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra 
    … …  
    566566 
    567567      !                             !* masked diffusive flux coefficients 
    568       ahu_bbl_0(:,:) = rn_ahtbbl * e2u(:,:) * e3u_bbl_0(:,:) / e1u(:,:) * umask(:,:,1) 
    569       ahv_bbl_0(:,:) = rn_ahtbbl * e1v(:,:) * e3v_bbl_0(:,:) / e2v(:,:) * vmask(:,:,1) 
     568      ahu_bbl_0(:,:) = rn_ahtbbl * e2_e1u(:,:) * e3u_bbl_0(:,:) * umask(:,:,1) 
     569      ahv_bbl_0(:,:) = rn_ahtbbl * e1_e2v(:,:) * e3v_bbl_0(:,:) * vmask(:,:,1) 
    570570 
    571571 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90

    r5102 r5836  
    66   !! History :  OPA  ! 1991-03  (O. Marti, G. Madec)  Original code 
    77   !!                 ! 1992-06  (M. Imbard)  doctor norme 
    8    !!                 ! 1996-01  (G. Madec)  statement function for e3 
    9    !!                 ! 1997-05  (G. Madec)  macro-tasked on jk-slab 
    108   !!                 ! 1998-07  (M. Imbard, G. Madec) ORCA version 
    11    !!            7.0  ! 2001-02  (M. Imbard)  cofdis, Original code 
     9   !!            7.0  ! 2001-02  (M. Imbard)  add distance to coast, Original code 
    1210   !!            8.1  ! 2001-02  (G. Madec, E. Durand)  cleaning 
    1311   !!  NEMO      1.0  ! 2002-08  (G. Madec, E. Durand)  free form + modules 
    … …  
    1513   !!            3.3  ! 2010-06  (C. Ethe, G. Madec) merge TRA-TRC  
    1614   !!            3.4  ! 2011-04  (G. Madec, C. Ethe) Merge of dtatem and dtasal + suppression of CPP keys 
     15   !!            3.6  ! 2015-06  (T. Graham)  read restoring coefficient in a file 
     16   !!            3.7  ! 2015-10  (G. Madec)  remove useless trends arrays 
    1717   !!---------------------------------------------------------------------- 
    1818 
    … …  
    4242 
    4343   PUBLIC   tra_dmp      ! routine called by step.F90 
    44    PUBLIC   tra_dmp_init ! routine called by opa.F90 
    45  
    46    !                               !!* Namelist namtra_dmp : T & S newtonian damping * 
    47    ! nn_zdmp and cn_resto are public as they are used by C1D/dyndmp.F90 
    48    LOGICAL , PUBLIC ::   ln_tradmp   !: internal damping flag 
    49    INTEGER , PUBLIC ::   nn_zdmp     ! = 0/1/2 flag for damping in the mixed layer 
    50    CHARACTER(LEN=200) , PUBLIC :: cn_resto      ! name of netcdf file containing restoration coefficient field 
     44   PUBLIC   tra_dmp_init ! routine called by nemogcm.F90 
     45 
     46   !                                           !!* Namelist namtra_dmp : T & S newtonian damping * 
     47   LOGICAL            , PUBLIC ::   ln_tradmp   !: internal damping flag 
     48   INTEGER            , PUBLIC ::   nn_zdmp     !: = 0/1/2 flag for damping in the mixed layer 
     49   CHARACTER(LEN=200) , PUBLIC ::   cn_resto    !: name of netcdf file containing restoration coefficient field 
    5150   ! 
    52  
    53  
    54    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   strdmp   !: damping salinity trend (psu/s) 
    55    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ttrdmp   !: damping temperature trend (Celcius/s) 
    5651   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   resto    !: restoring coeff. on T and S (s-1) 
    5752 
    … …  
    7065      !!                ***  FUNCTION tra_dmp_alloc  *** 
    7166      !!---------------------------------------------------------------------- 
    72       ALLOCATE( strdmp(jpi,jpj,jpk) , ttrdmp(jpi,jpj,jpk), resto(jpi,jpj,jpk), STAT= tra_dmp_alloc ) 
     67      ALLOCATE( resto(jpi,jpj,jpk), STAT= tra_dmp_alloc ) 
    7368      ! 
    7469      IF( lk_mpp            )   CALL mpp_sum ( tra_dmp_alloc ) 
    … …  
    9691      !! ** Action  : - (ta,sa)   tracer trends updated with the damping trend 
    9792      !!---------------------------------------------------------------------- 
    98       ! 
    9993      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    100       !! 
    101       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    102       REAL(wp) ::   zta, zsa             ! local scalars 
    103       REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zts_dta  
    104       !!---------------------------------------------------------------------- 
    105       ! 
    106       IF( nn_timing == 1 )  CALL timing_start( 'tra_dmp') 
    107       ! 
    108       CALL wrk_alloc( jpi, jpj, jpk, jpts,  zts_dta ) 
    109       ! 
    110       !                           !==   input T-S data at kt   ==! 
     94      ! 
     95      INTEGER ::   ji, jj, jk, jn   ! dummy loop indices 
     96      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zts_dta, ztrdts 
     97      !!---------------------------------------------------------------------- 
     98      ! 
     99      IF( nn_timing == 1 )   CALL timing_start('tra_dmp') 
     100      ! 
     101      CALL wrk_alloc( jpi,jpj,jpk,jpts,   zts_dta ) 
     102      ! 
     103      IF( l_trdtra )   THEN                    !* Save ta and sa trends 
     104         CALL wrk_alloc( jpi,jpj,jpk,jpts,   ztrdts )  
     105         ztrdts(:,:,:,:) = tsa(:,:,:,:)  
     106      ENDIF 
     107      !                           !==  input T-S data at kt  ==! 
    111108      CALL dta_tsd( kt, zts_dta )            ! read and interpolates T-S data at kt 
    112109      ! 
    113       SELECT CASE ( nn_zdmp )     !==    type of damping   ==! 
    114       ! 
    115       CASE( 0 )                   !==  newtonian damping throughout the water column  ==! 
    116          DO jk = 1, jpkm1 
    117             DO jj = 2, jpjm1 
    118                DO ji = fs_2, fs_jpim1   ! vector opt. 
    119                   zta = resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) 
    120                   zsa = resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 
    121                   tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta 
    122                   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa 
    123                   strdmp(ji,jj,jk) = zsa           ! save the trend (used in asmtrj) 
    124                   ttrdmp(ji,jj,jk) = zta       
     110      SELECT CASE ( nn_zdmp )     !==  type of damping  ==! 
     111      ! 
     112      CASE( 0 )                        !*  newtonian damping throughout the water column  *! 
     113         DO jn = 1, jpts 
     114            DO jk = 1, jpkm1 
     115               DO jj = 2, jpjm1 
     116                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     117                     tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - tsb(ji,jj,jk,jn) ) 
     118                  END DO 
    125119               END DO 
    126120            END DO 
    127121         END DO 
    128122         ! 
    129       CASE ( 1 )                  !==  no damping in the turbocline (avt > 5 cm2/s)  ==! 
     123      CASE ( 1 )                       !*  no damping in the turbocline (avt > 5 cm2/s)  *! 
    130124         DO jk = 1, jpkm1 
    131125            DO jj = 2, jpjm1 
    132126               DO ji = fs_2, fs_jpim1   ! vector opt. 
    133127                  IF( avt(ji,jj,jk) <= 5.e-4_wp ) THEN 
    134                      zta = resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) 
    135                      zsa = resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 
    136                   ELSE 
    137                      zta = 0._wp 
    138                      zsa = 0._wp   
     128                     tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   & 
     129                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) 
     130                     tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   & 
     131                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 
    139132                  ENDIF 
    140                   tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta 
    141                   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa 
    142                   strdmp(ji,jj,jk) = zsa           ! save the salinity trend (used in asmtrj) 
    143                   ttrdmp(ji,jj,jk) = zta 
    144133               END DO 
    145134            END DO 
    146135         END DO 
    147136         ! 
    148       CASE ( 2 )                  !==  no damping in the mixed layer   ==! 
     137      CASE ( 2 )                       !*  no damping in the mixed layer   *! 
    149138         DO jk = 1, jpkm1 
    150139            DO jj = 2, jpjm1 
    151140               DO ji = fs_2, fs_jpim1   ! vector opt. 
    152141                  IF( fsdept(ji,jj,jk) >= hmlp (ji,jj) ) THEN 
    153                      zta = resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) 
    154                      zsa = resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 
    155                   ELSE 
    156                      zta = 0._wp 
    157                      zsa = 0._wp   
     142                     tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   & 
     143                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) 
     144                     tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   & 
     145                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 
    158146                  ENDIF 
    159                   tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta 
    160                   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa 
    161                   strdmp(ji,jj,jk) = zsa           ! save the salinity trend (used in asmtrj) 
    162                   ttrdmp(ji,jj,jk) = zta 
    163147               END DO 
    164148            END DO 
    … …  
    168152      ! 
    169153      IF( l_trdtra )   THEN       ! trend diagnostic 
    170          CALL trd_tra( kt, 'TRA', jp_tem, jptra_dmp, ttrdmp ) 
    171          CALL trd_tra( kt, 'TRA', jp_sal, jptra_dmp, strdmp ) 
     154         ztrdts(:,:,:,:) = tsa(:,:,:,:) - ztrdts(:,:,:,:) 
     155         CALL trd_tra( kt, 'TRA', jp_tem, jptra_dmp, ztrdts(:,:,:,jp_tem) ) 
     156         CALL trd_tra( kt, 'TRA', jp_sal, jptra_dmp, ztrdts(:,:,:,jp_sal) ) 
     157         CALL wrk_dealloc( jpi,jpj,jpk,jpts,   ztrdts )  
    172158      ENDIF 
    173159      !                           ! Control print 
    … …  
    175161         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    176162      ! 
    177       CALL wrk_dealloc( jpi, jpj, jpk, jpts,  zts_dta ) 
    178       ! 
    179       IF( nn_timing == 1 )  CALL timing_stop( 'tra_dmp') 
     163      CALL wrk_dealloc( jpi,jpj,jpk,jpts,   zts_dta ) 
     164      ! 
     165      IF( nn_timing == 1 )   CALL timing_stop('tra_dmp') 
    180166      ! 
    181167   END SUBROUTINE tra_dmp 
    … …  
    190176      !! ** Method  :   read the namtra_dmp namelist and check the parameters 
    191177      !!---------------------------------------------------------------------- 
     178      INTEGER ::   ios, imask   ! local integers  
     179      !! 
    192180      NAMELIST/namtra_dmp/ ln_tradmp, nn_zdmp, cn_resto 
    193       INTEGER ::  ios         ! Local integer for output status of namelist read 
    194       INTEGER :: imask        ! File handle  
    195       !! 
    196181      !!---------------------------------------------------------------------- 
    197182      ! 
    … …  
    204189902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_dmp in configuration namelist', lwp ) 
    205190      IF(lwm) WRITE ( numond, namtra_dmp ) 
    206  
    207       IF(lwp) THEN                 !Namelist print 
     191      ! 
     192      IF(lwp) THEN                  ! Namelist print 
    208193         WRITE(numout,*) 
    209194         WRITE(numout,*) 'tra_dmp_init : T and S newtonian relaxation' 
    210          WRITE(numout,*) '~~~~~~~' 
     195         WRITE(numout,*) '~~~~~~~~~~~' 
    211196         WRITE(numout,*) '   Namelist namtra_dmp : set relaxation parameters' 
    212197         WRITE(numout,*) '      Apply relaxation   or not       ln_tradmp = ', ln_tradmp 
    … …  
    215200         WRITE(numout,*) 
    216201      ENDIF 
    217  
     202      ! 
    218203      IF( ln_tradmp) THEN 
    219          ! 
    220          !Allocate arrays 
     204         !                          ! Allocate arrays 
    221205         IF( tra_dmp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_dmp_init: unable to allocate arrays' ) 
    222  
    223          !Check values of nn_zdmp 
    224          SELECT CASE (nn_zdmp) 
    225          CASE ( 0 )  ; IF(lwp) WRITE(numout,*) '   tracer damping as specified by mask' 
    226          CASE ( 1 )  ; IF(lwp) WRITE(numout,*) '   no tracer damping in the turbocline' 
    227          CASE ( 2 )  ; IF(lwp) WRITE(numout,*) '   no tracer damping in the mixed layer' 
     206         ! 
     207         SELECT CASE (nn_zdmp)      ! Check values of nn_zdmp 
     208         CASE ( 0 )   ;   IF(lwp) WRITE(numout,*) '   tracer damping as specified by mask' 
     209         CASE ( 1 )   ;   IF(lwp) WRITE(numout,*) '   no tracer damping in the mixing layer (kz > 5 cm2/s)' 
     210         CASE ( 2 )   ;   IF(lwp) WRITE(numout,*) '   no tracer damping in the mixed  layer' 
     211         CASE DEFAULT 
     212            CALL ctl_stop('tra_dmp_init : wrong value of nn_zdmp') 
    228213         END SELECT 
    229  
    230          !TG: Initialisation of dtatsd - Would it be better to have dmpdta routine 
    231          !so can damp to something other than intitial conditions files? 
     214         ! 
     215         !!TG: Initialisation of dtatsd - Would it be better to have dmpdta routine 
     216         !    so can damp to something other than intitial conditions files? 
     217         !!gm: In principle yes. Nevertheless, we can't anticipate demands that have never been formulated. 
    232218         IF( .NOT.ln_tsd_tradmp ) THEN 
    233             CALL ctl_warn( 'tra_dmp_init: read T-S data not initialized, we force ln_tsd_tradmp=T' ) 
     219            IF(lwp) WRITE(numout,*) 
     220            IF(lwp) WRITE(numout, *)  '   read T-S data not initialized, we force ln_tsd_tradmp=T' 
    234221            CALL dta_tsd_init( ld_tradmp=ln_tradmp )        ! forces the initialisation of T-S data 
    235222         ENDIF 
    236  
    237          !initialise arrays - Are these actually used anywhere else? 
    238          strdmp(:,:,:) = 0._wp 
    239          ttrdmp(:,:,:) = 0._wp 
    240  
    241          !Read in mask from file 
     223         !                          ! Read in mask from file 
    242224         CALL iom_open ( cn_resto, imask) 
    243          CALL iom_get  ( imask, jpdom_autoglo, 'resto', resto) 
     225         CALL iom_get  ( imask, jpdom_autoglo, 'resto', resto ) 
    244226         CALL iom_close( imask ) 
    245        ENDIF 
    246  
     227      ENDIF 
     228      ! 
    247229   END SUBROUTINE tra_dmp_init 
    248230 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90

    r5120 r5836  
    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   ! 
    41    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  
     42   INTEGER ::   nldf = 0   ! type of lateral diffusion used defined from ln_traldf_... (namlist logicals) 
     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. 
    90          ! 
    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' ) 
     72      ! 
     73      SELECT CASE ( nldf )                     !* compute lateral mixing trend and add it to the general trend 
     74      ! 
     75      CASE ( np_lap   )                                  ! laplacian: iso-level operator 
     76         CALL tra_ldf_lap  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb,      tsa, jpts,  1   ) 
     77      CASE ( np_lap_i )                                  ! laplacian: standard iso-neutral operator (Madec) 
     78         CALL tra_ldf_iso  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts,  1   ) 
     79      CASE ( np_lap_it )                                 ! laplacian: triad iso-neutral operator (griffies) 
     80         CALL tra_ldf_triad( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts,  1   ) 
     81      CASE ( np_blp , np_blp_i , np_blp_it )             ! bilaplacian: iso-level & iso-neutral operators 
     82         CALL tra_ldf_blp  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb      , tsa, jpts, nldf ) 
    11183      END SELECT 
    11284 
    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 
     85      IF( l_trdtra )   THEN                    !* save the horizontal diffusive trends for further diagnostics 
    11986         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    12087         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
    12188         CALL trd_tra( kt, 'TRA', jp_tem, jptra_ldf, ztrdt ) 
    12289         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) 
     90         CALL wrk_dealloc( jpi,jpj,jpk,  ztrdt, ztrds )  
     91      ENDIF 
     92      !                                        !* print mean trends (used for debugging) 
    12693      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf  - Ta: ', mask1=tmask,               & 
    12794         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    12895      ! 
    129       IF( nn_timing == 1 )  CALL timing_stop('tra_ldf') 
     96      IF( nn_timing == 1 )   CALL timing_stop('tra_ldf') 
    13097      ! 
    13198   END SUBROUTINE tra_ldf 
    … …  
    139106      !! 
    140107      !! ** 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 
     108      !!---------------------------------------------------------------------- 
     109      INTEGER ::   ioptio, ierr   ! temporary integers  
     110      !!---------------------------------------------------------------------- 
     111      ! 
     112      IF(lwp) THEN                     ! Namelist print 
    154113         WRITE(numout,*) 
    155114         WRITE(numout,*) 'tra_ldf_init : lateral tracer diffusive operator' 
    … …  
    159118         WRITE(numout,*) 
    160119      ENDIF 
    161  
    162       !                               ! control the input 
     120      !                                   ! use of lateral operator or not 
     121      nldf   = np_ERROR 
    163122      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 
    168       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 
    176       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) 
     123      IF( ln_traldf_lap )   ioptio = ioptio + 1 
     124      IF( ln_traldf_blp )   ioptio = ioptio + 1 
     125      IF( ioptio >  1   )   CALL ctl_stop( 'tra_ldf_init: use ONE or NONE of the 2 lap/bilap operator type on tracer' ) 
     126      IF( ioptio == 0   )   nldf = np_no_ldf     ! No lateral diffusion 
     127      ! 
     128      IF( nldf /= np_no_ldf ) THEN        ! direction ==>> type of operator   
     129         ioptio = 0 
     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 
     136         ierr = 0 
     137         IF( ln_traldf_lap ) THEN         ! laplacian operator 
     138            IF ( ln_zco ) THEN               ! z-coordinate 
     139               IF ( ln_traldf_lev   )   nldf = np_lap     ! iso-level = horizontal (no rotation) 
     140               IF ( ln_traldf_hor   )   nldf = np_lap     ! iso-level = horizontal (no rotation) 
     141               IF ( ln_traldf_iso   )   nldf = np_lap_i   ! iso-neutral: standard  (   rotation) 
     142               IF ( ln_traldf_triad )   nldf = np_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 = np_lap     ! horizontal             (no rotation) 
     147               IF ( ln_traldf_iso   )   nldf = np_lap_i   ! iso-neutral: standard     (rotation) 
     148               IF ( ln_traldf_triad )   nldf = np_lap_it  ! iso-neutral: triad        (rotation) 
     149            ENDIF 
     150            IF ( ln_sco ) THEN               ! s-coordinate 
     151               IF ( ln_traldf_lev   )   nldf = np_lap     ! iso-level              (no rotation) 
     152               IF ( ln_traldf_hor   )   nldf = np_lap_i   ! horizontal             (   rotation) 
     153               IF ( ln_traldf_iso   )   nldf = np_lap_i   ! iso-neutral: standard  (   rotation) 
     154               IF ( ln_traldf_triad )   nldf = np_lap_it  ! iso-neutral: triad     (   rotation) 
     155            ENDIF 
    182156         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) 
     157         ! 
     158         IF( ln_traldf_blp ) THEN         ! bilaplacian operator 
     159            IF ( ln_zco ) THEN               ! z-coordinate 
     160               IF ( ln_traldf_lev   )   nldf = np_blp     ! iso-level = horizontal (no rotation) 
     161               IF ( ln_traldf_hor   )   nldf = np_blp     ! iso-level = horizontal (no rotation) 
     162               IF ( ln_traldf_iso   )   nldf = np_blp_i   ! iso-neutral: standard  (   rotation) 
     163               IF ( ln_traldf_triad )   nldf = np_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 = np_blp     ! horizontal             (no rotation) 
     168               IF ( ln_traldf_iso   )   nldf = np_blp_i   ! iso-neutral: standard  (   rotation) 
     169               IF ( ln_traldf_triad )   nldf = np_blp_it  ! iso-neutral: triad     (   rotation) 
     170            ENDIF 
     171            IF ( ln_sco ) THEN               ! s-coordinate 
     172               IF ( ln_traldf_lev   )   nldf = np_blp     ! iso-level              (no rotation) 
     173               IF ( ln_traldf_hor   )   nldf = np_blp_it  ! horizontal             (   rotation) 
     174               IF ( ln_traldf_iso   )   nldf = np_blp_i   ! iso-neutral: standard  (   rotation) 
     175               IF ( ln_traldf_triad )   nldf = np_blp_it  ! iso-neutral: triad     (   rotation) 
     176            ENDIF 
    187177         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' ) 
     178      ENDIF 
     179      ! 
    214180      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  
     181      IF( ln_ldfeiv .AND. .NOT.( ln_traldf_iso .OR. ln_traldf_triad ) )                                    & 
     182           &            CALL ctl_stop( '          eddy induced velocity on tracers requires isopycnal',    & 
     183           &                                                                    ' laplacian diffusion' ) 
     184      IF(  nldf == np_lap_i .OR. nldf == np_lap_it .OR. & 
     185         & nldf == np_blp_i .OR. nldf == np_blp_it  )   l_ldfslp = .TRUE.    ! slope of neutral surfaces required  
     186      ! 
    229187      IF(lwp) THEN 
    230188         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 
     189         IF( nldf == np_no_ldf )   WRITE(numout,*) '          NO lateral diffusion' 
     190         IF( nldf == np_lap    )   WRITE(numout,*) '          laplacian iso-level operator' 
     191         IF( nldf == np_lap_i  )   WRITE(numout,*) '          Rotated laplacian operator (standard)' 
     192         IF( nldf == np_lap_it )   WRITE(numout,*) '          Rotated laplacian operator (triad)' 
     193         IF( nldf == np_blp    )   WRITE(numout,*) '          bilaplacian iso-level operator' 
     194         IF( nldf == np_blp_i  )   WRITE(numout,*) '          Rotated bilaplacian operator (standard)' 
     195         IF( nldf == np_blp_it )   WRITE(numout,*) '          Rotated bilaplacian operator (triad)' 
     196      ENDIF 
    242197      ! 
    243198   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 
    361199 
    362200   !!====================================================================== 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90

    r5149 r5836  
    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 ) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk) 
    203                   zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * re1v_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 / ( e12t(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 / ( e12t(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   !!============================================================================== 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90

    r5147 r5836  
    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) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk) 
    99                   zabe2 = fsahtv(ji,jj,jk) * vmask(ji,jj,jk) * re1v_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) * re2u_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) * re1v_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) * re2u_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) * re1v_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 / ( e12t(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 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90

    r5656 r5836  
    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) 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

    r5407 r5836  
    189189                  CALL fld_read( kt, 1, sf_chl )                         ! Read Chl data and provides it at the current time step 
    190190                  !          
    191 !CDIR COLLAPSE 
    192 !CDIR NOVERRCHK 
    193191                  DO jj = 1, jpj                                         ! Separation in R-G-B depending of the surface Chl 
    194 !CDIR NOVERRCHK 
    195192                     DO ji = 1, jpi 
    196193                        zchl = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) 
    … …  
    217214               ! 
    218215               DO jk = 2, nksr+1 
    219 !CDIR NOVERRCHK 
    220216                  DO jj = 1, jpj 
    221 !CDIR NOVERRCHK    
    222217                     DO ji = 1, jpi 
    223218                        zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r     ) 
    … …  
    495490                
    496491                  DO jk = 2, nksr+1 
    497 !CDIR NOVERRCHK 
    498492                     DO jj = 1, jpj 
    499 !CDIR NOVERRCHK    
    500493                        DO ji = 1, jpi 
    501494                           zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * xsi0r     ) 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90

    r5385 r5836  
    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) 
    … …  
    8083      CASE ( 0 )    ;    CALL tra_zdf_exp( kt, nit000, 'TRA', r2dtra, nn_zdfexp, tsb, tsa, jpts )  !   explicit scheme  
    8184      CASE ( 1 )    ;    CALL tra_zdf_imp( kt, nit000, 'TRA', r2dtra,            tsb, tsa, jpts )  !   implicit scheme  
    82       CASE ( -1 )                                       ! esopa: test all possibility with control print 
    83          CALL tra_zdf_exp( kt, nit000, 'TRA', r2dtra, nn_zdfexp, tsb, tsa, jpts ) 
    84          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' zdf0 - Ta: ', mask1=tmask,               & 
    85          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    86          CALL tra_zdf_imp( kt, nit000, 'TRA', r2dtra,            tsb, tsa, jpts )  
    87          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' zdf1 - Ta: ', mask1=tmask,               & 
    88          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    8985      END SELECT 
     86!!gm WHY here !   and I don't like that ! 
    9087      ! DRAKKAR SSS control { 
    9188      ! JMM avoid negative salinities near river outlet ! Ugly fix 
    9289      ! JMM : restore negative salinities to small salinities: 
    9390      WHERE ( tsa(:,:,:,jp_sal) < 0._wp )   tsa(:,:,:,jp_sal) = 0.1_wp 
     91!!gm 
    9492 
    9593      IF( l_trdtra )   THEN                      ! save the vertical diffusive trends for further diagnostics 
    … …  
    9896            ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal) - tsb(:,:,jk,jp_sal) ) / r2dtra(jk) ) - ztrds(:,:,jk) 
    9997         END DO 
     98!!gm this should be moved in trdtra.F90 and done on all trends 
    10099         CALL lbc_lnk( ztrdt, 'T', 1. ) 
    101100         CALL lbc_lnk( ztrds, 'T', 1. ) 
     101!!gm 
    102102         CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt ) 
    103103         CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds ) 
    … …  
    123123      !!      nzdf = 0   explicit (time-splitting) scheme (ln_zdfexp=T) 
    124124      !!           = 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. 
     125      !!      NB: rotation of lateral mixing operator or TKE & GLS schemes, 
     126      !!          an implicit scheme is required. 
    127127      !!---------------------------------------------------------------------- 
    128128      USE zdftke 
    129129      USE zdfgls 
    130       USE zdfkpp 
    131130      !!---------------------------------------------------------------------- 
    132131 
    … …  
    137136 
    138137      ! 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 
     138      IF( lk_zdftke .OR. lk_zdfgls   )   nzdf = 1   ! TKE, or GLS physics 
     139      IF( ln_traldf_iso              )   nzdf = 1   ! iso-neutral lateral physics 
     140      IF( ln_traldf_hor .AND. ln_sco )   nzdf = 1   ! horizontal lateral physics in s-coordinate 
    142141      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.' ) 
    144  
    145       ! Test: esopa 
    146       IF( lk_esopa )    nzdf = -1                      ! All schemes used 
     142            &                         ' GLS or TKE scheme, the implicit scheme is required, set ln_zdfexp = .false.' ) 
    147143 
    148144      IF(lwp) THEN 
    … …  
    150146         WRITE(numout,*) 'tra_zdf_init : vertical tracer physics scheme' 
    151147         WRITE(numout,*) '~~~~~~~~~~~' 
    152          IF( nzdf == -1 )   WRITE(numout,*) '              ESOPA test All scheme used' 
    153148         IF( nzdf ==  0 )   WRITE(numout,*) '              Explicit time-splitting scheme' 
    154149         IF( nzdf ==  1 )   WRITE(numout,*) '              Implicit (euler backward) scheme' 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90

    r5120 r5836  
    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_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) 
    … …  
    7674      !! ** Action  : - pta  becomes the after tracer 
    7775      !!--------------------------------------------------------------------- 
    78       USE oce     , ONLY:   zwd => ua       , zws => va         ! (ua,va) used as 3D workspace 
    79       ! 
    8076      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index 
    8177      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index 
    … …  
    8884      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
    8985      REAL(wp) ::  zrhs, ze3tb, ze3tn, ze3ta   ! local scalars 
    90       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwt 
     86      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwt, zwd, zws 
    9187      !!--------------------------------------------------------------------- 
    9288      ! 
    9389      IF( nn_timing == 1 )  CALL timing_start('tra_zdf_imp') 
    9490      ! 
    95       CALL wrk_alloc( jpi, jpj, jpk, zwi, zwt )  
     91      CALL wrk_alloc( jpi,jpj,jpk,   zwi, zwt, zwd, zws )  
    9692      ! 
    9793      IF( kt == kit000 )  THEN 
    … …  
    120116            ELSE                                            ;   zwt(:,:,2:jpk) = fsavs(:,:,2:jpk) 
    121117            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)        
     118            zwt(:,:,1) = 0._wp 
     119            ! 
     120            IF( l_ldfslp ) THEN            ! isoneutral diffusion: add the contribution  
     121               IF( ln_traldf_msc  ) THEN     ! MSC iso-neutral operator  
     122                  DO jk = 2, jpkm1 
     123                     DO jj = 2, jpjm1 
     124                        DO ji = fs_2, fs_jpim1   ! vector opt. 
     125                           zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk)   
     126                        END DO 
    135127                     END DO 
    136128                  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)  ) 
     129               ELSE                          ! standard or triad iso-neutral operator 
     130                  DO jk = 2, jpkm1 
     131                     DO jj = 2, jpjm1 
     132                        DO ji = fs_2, fs_jpim1   ! vector opt. 
     133                           zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk) 
     134                        END DO 
    145135                     END DO 
    146136                  END DO 
    147                END DO 
     137               ENDIF 
    148138            ENDIF 
    149 #endif 
     139            ! 
    150140            ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked) 
    151141            DO jk = 1, jpkm1 
    … …  
    202192               ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,1) 
    203193               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) 
     194               pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) + p2dt(1) * ze3tn * pta(ji,jj,1,jn) 
    206195            END DO 
    207196         END DO 
    … …  
    235224      !                                               ! ================= ! 
    236225      ! 
    237       CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwt )  
     226      CALL wrk_dealloc( jpi,jpj,jpk,   zwi, zwt, zwd, zws )  
    238227      ! 
    239228      IF( nn_timing == 1 )  CALL timing_stop('tra_zdf_imp') 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90

    r5120 r5836  
    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.