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 5972 for branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90 – NEMO

Ignore:
Timestamp:
2015-12-02T09:52:20+01:00 (8 years ago)
Author:
timgraham
Message:

Upgraded to head of trunk (r5936)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90

    r5967 r5972  
    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  
    19    USE dynspg_oce     ! choice/control of key cpp for surface pressure gradient 
    2020   USE diaptr         ! poleward transport diagnostics 
    2121   ! 
     
    3838#  include "vectopt_loop_substitute.h90" 
    3939   !!---------------------------------------------------------------------- 
    40    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     40   !! NEMO/OPA 3.7 , NEMO Consortium (2015) 
    4141   !! $Id$ 
    4242   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    4444CONTAINS 
    4545 
    46    SUBROUTINE tra_adv_ubs ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      & 
    47       &                                       ptb, ptn, pta, kjpt ) 
     46   SUBROUTINE tra_adv_ubs( kt, kit000, cdtype, p2dt, pun, pvn, pwn,          & 
     47      &                                                ptb, ptn, pta, kjpt, kn_ubs_v ) 
    4848      !!---------------------------------------------------------------------- 
    4949      !!                  ***  ROUTINE tra_adv_ubs  *** 
     
    5252      !!      and add it to the general trend of passive tracer equations. 
    5353      !! 
    54       !! ** Method  :   The upstream biased scheme (UBS) is based on a 3rd order 
     54      !! ** Method  :   The 3rd order Upstream Biased Scheme (UBS) is based on an 
    5555      !!      upstream-biased parabolic interpolation (Shchepetkin and McWilliams 2005) 
    5656      !!      It is only used in the horizontal direction. 
     
    6161      !!      where zltu is the second derivative of the before temperature field: 
    6262      !!          zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ] 
    63       !!      This results in a dissipatively dominant (i.e. hyper-diffusive)  
     63      !!        This results in a dissipatively dominant (i.e. hyper-diffusive)  
    6464      !!      truncation error. The overall performance of the advection scheme  
    6565      !!      is similar to that reported in (Farrow and Stevens, 1995).  
    66       !!      For stability reasons, the first term of the fluxes which corresponds 
     66      !!        For stability reasons, the first term of the fluxes which corresponds 
    6767      !!      to a second order centered scheme is evaluated using the now velocity  
    6868      !!      (centered in time) while the second term which is the diffusive part  
    6969      !!      of the scheme, is evaluated using the before velocity (forward in time).  
    7070      !!      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. 
     71      !!                On the vertical, the advection is evaluated using a FCT scheme, 
     72      !!      as the UBS have been found to be too diffusive.  
     73!!gm  !!                kn_ubs_v argument (not coded for the moment) 
     74      !!      controles whether the FCT is based on a 2nd order centrered scheme (kn_ubs_v=2)  
     75      !!      or on a 4th order compact scheme (kn_ubs_v=4). 
    7376      !! 
    7477      !! ** Action : - update (pta) with the now advective tracer trends 
     
    8184      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator) 
    8285      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers 
     86      INTEGER                              , INTENT(in   ) ::   kn_ubs_v        ! number of tracers 
    8387      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step 
    8488      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean transport components 
     
    9599      IF( nn_timing == 1 )  CALL timing_start('tra_adv_ubs') 
    96100      ! 
    97       CALL wrk_alloc( jpi, jpj, jpk, ztu, ztv, zltu, zltv, zti, ztw ) 
     101      CALL wrk_alloc( jpi,jpj,jpk,  ztu, ztv, zltu, zltv, zti, ztw ) 
    98102      ! 
    99103      IF( kt == kit000 )  THEN 
     
    106110      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 
    107111      ! 
     112      zltu(:,:,jpk) = 0._wp   ;   zltv(:,:,jpk) = 0._wp     ! Bottom value : set to zero one for all 
     113      ztw (:,:,jpk) = 0._wp   ;   zti (:,:,jpk) = 0._wp 
     114      IF( lk_vvl )   ztw(:,:, 1 ) = 0._wp                   ! surface value: set to zero only in vvl case 
     115      ! 
    108116      !                                                          ! =========== 
    109117      DO jn = 1, kjpt                                            ! tracer loop 
    110118         !                                                       ! =========== 
    111          ! 1. Bottom value : flux set to zero 
    112          ! ---------------------------------- 
    113          zltu(:,:,jpk) = 0.e0       ;      zltv(:,:,jpk) = 0.e0 
    114119         !                                               
    115          DO jk = 1, jpkm1                                 ! Horizontal slab 
    116             !                                    
    117             !  Laplacian 
    118             DO jj = 1, jpjm1            ! First derivative (gradient) 
     120         DO jk = 1, jpkm1        !==  horizontal laplacian of before tracer ==! 
     121            DO jj = 1, jpjm1              ! First derivative (masked gradient) 
    119122               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) 
     123                  zeeu = e2_e1u(ji,jj) * fse3u(ji,jj,jk) * umask(ji,jj,jk) 
     124                  zeev = e1_e2v(ji,jj) * fse3v(ji,jj,jk) * vmask(ji,jj,jk) 
    122125                  ztu(ji,jj,jk) = zeeu * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) 
    123126                  ztv(ji,jj,jk) = zeev * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) 
    124127               END DO 
    125128            END DO 
    126             DO jj = 2, jpjm1            ! Second derivative (divergence) 
     129            DO jj = 2, jpjm1              ! Second derivative (divergence) 
    127130               DO ji = fs_2, fs_jpim1   ! vector opt. 
    128                   zcoef = 1. / ( 6. * fse3t(ji,jj,jk) ) 
     131                  zcoef = 1._wp / ( 6._wp * fse3t(ji,jj,jk) ) 
    129132                  zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef 
    130133                  zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef 
     
    132135            END DO 
    133136            !                                     
    134          END DO                                           ! End of slab          
     137         END DO          
    135138         CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn) 
    136  
    137139         !     
    138          !  Horizontal advective fluxes                
    139          DO jk = 1, jpkm1                                 ! Horizontal slab 
     140         DO jk = 1, jpkm1        !==  Horizontal advective fluxes  ==!     (UBS) 
    140141            DO jj = 1, jpjm1 
    141142               DO ji = 1, fs_jpim1   ! vector opt. 
    142                   ! upstream transport (x2) 
    143                   zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) 
     143                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )      ! upstream transport (x2) 
    144144                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) 
    145145                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) 
    146146                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) 
    147                   ! 2nd order centered advective fluxes (x2) 
     147                  !                                                  ! 2nd order centered advective fluxes (x2) 
    148148                  zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn) ) 
    149149                  zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) ) 
    150                   ! UBS advective fluxes 
     150                  !                                                  ! UBS advective fluxes 
    151151                  ztu(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) ) 
    152152                  ztv(ji,jj,jk) = 0.5 * ( zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) ) 
    153153               END DO 
    154154            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 
     155         END DO          
     156         ! 
     157         zltu(:,:,:) = pta(:,:,:,jn)      ! store the initial trends before its update 
     158         ! 
     159         DO jk = 1, jpkm1        !==  add the horizontal advective trend  ==! 
    160160            DO jj = 2, jpjm1 
    161161               DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    166166            END DO 
    167167            !                                              
    168          END DO                                           !   End of slab 
    169  
    170          ! Horizontal trend used in tra_adv_ztvd subroutine 
    171          zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:) 
    172  
     168         END DO 
     169         ! 
     170         zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:)    ! Horizontal advective trend used in vertical 2nd order FCT case 
     171         !                                            ! and/or in trend diagnostic (l_trd=T)  
    173172         !                 
    174173         IF( l_trd ) THEN                  ! trend diagnostics 
     
    181180            IF( jn == jp_sal )  str_adv(:) = ptr_sj( ztv(:,:,:) ) 
    182181         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 
     182         ! 
     183         !                       !== vertical advective trend  ==! 
     184         ! 
     185         SELECT CASE( kn_ubs_v )       ! select the vertical advection scheme 
     186         ! 
     187         CASE(  2  )                   ! 2nd order FCT  
     188            !          
     189            IF( l_trd )   zltv(:,:,:) = pta(:,:,:,jn)          ! store pta if trend diag. 
     190            ! 
     191            !                          !*  upstream advection with initial mass fluxes & intermediate update  ==! 
     192            DO jk = 2, jpkm1                 ! Interior value (w-masked) 
     193               DO jj = 1, jpj 
     194                  DO ji = 1, jpi 
     195                     zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 
     196                     zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 
     197                     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) 
     198                  END DO 
     199               END DO 
     200            END DO  
     201            IF(.NOT.lk_vvl ) THEN            ! top ocean value (only in linear free surface as ztw has been w-masked) 
     202               IF( ln_isfcav ) THEN                ! top of the ice-shelf cavities and at the ocean surface 
     203                  DO jj = 1, jpj 
     204                     DO ji = 1, jpi 
     205                        ztw(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
     206                     END DO 
     207                  END DO    
     208               ELSE                                ! no cavities: only at the ocean surface 
     209                  ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 
     210               ENDIF 
     211            ENDIF 
     212            ! 
     213            DO jk = 1, jpkm1           !* trend and after field with monotonic scheme 
     214               z2dtt = p2dt(jk) 
     215               DO jj = 2, jpjm1 
     216                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     217                     ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 
     218                     pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn) +  ztak  
     219                     zti(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk) 
     220                  END DO 
     221               END DO 
     222            END DO 
     223            CALL lbc_lnk( zti, 'T', 1. )      ! Lateral boundary conditions on zti, zsi   (unchanged sign) 
     224            ! 
     225            !                          !*  anti-diffusive flux : high order minus low order 
     226            DO jk = 2, jpkm1        ! Interior value  (w-masked) 
     227               DO jj = 1, jpj 
     228                  DO ji = 1, jpi 
     229                     ztw(ji,jj,jk) = (   0.5_wp * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )   & 
     230                        &              - ztw(ji,jj,jk)   ) * wmask(ji,jj,jk) 
     231                  END DO 
     232               END DO 
     233            END DO 
     234            !                                            ! top ocean value: high order == upstream  ==>>  zwz=0 
     235            IF(.NOT.lk_vvl )   ztw(:,:, 1 ) = 0._wp      ! only ocean surface as interior zwz values have been w-masked 
     236            ! 
     237            CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt )      !  monotonicity algorithm 
     238            ! 
     239         CASE(  4  )                               ! 4th order COMPACT 
     240            CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )         ! 4th order compact interpolation of T at w-point 
     241            DO jk = 2, jpkm1 
     242               DO jj = 2, jpjm1 
     243                  DO ji = fs_2, fs_jpim1 
     244                     ztw(ji,jj,jk) = pwn(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk) 
     245                  END DO 
     246               END DO 
     247            END DO 
     248            IF(.NOT.lk_vvl )   ztw(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn)     !!gm ISF & 4th COMPACT doesn't work 
     249            ! 
     250         END SELECT 
     251         ! 
     252         DO jk = 1, jpkm1        !  final trend with corrected fluxes 
    236253            DO jj = 2, jpjm1  
    237254               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 
     255                  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) ) 
     256               END DO 
     257            END DO 
     258         END DO 
     259         ! 
     260         IF( l_trd )  THEN       ! vertical advective trend diagnostics 
    249261            DO jk = 1, jpkm1                       ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w]) 
    250262               DO jj = 2, jpjm1 
    251263                  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 
     264                     zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk)                          & 
     265                        &           + ptn(ji,jj,jk,jn) * (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  )   & 
     266                        &                              / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 
    255267                  END DO 
    256268               END DO 
     
    261273      END DO 
    262274      ! 
    263       CALL wrk_dealloc( jpi, jpj, jpk, ztu, ztv, zltu, zltv, zti, ztw ) 
     275      CALL wrk_dealloc( jpi,jpj,jpk,  ztu, ztv, zltu, zltv, zti, ztw ) 
    264276      ! 
    265277      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_ubs') 
     
    294306      IF( nn_timing == 1 )  CALL timing_start('nonosc_z') 
    295307      ! 
    296       CALL wrk_alloc( jpi, jpj, jpk, zbetup, zbetdo ) 
     308      CALL wrk_alloc( jpi,jpj,jpk,  zbetup, zbetdo ) 
    297309      ! 
    298310      zbig  = 1.e+40_wp 
    299311      zrtrn = 1.e-15_wp 
    300312      zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp 
    301  
     313      ! 
    302314      ! Search local extrema 
    303315      ! -------------------- 
    304       ! large negative value (-zbig) inside land 
     316      !                    ! large negative value (-zbig) inside land 
    305317      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) ) 
    306318      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) ) 
    307       ! search maximum in neighbourhood 
    308       DO jk = 1, jpkm1 
     319      ! 
     320      DO jk = 1, jpkm1     ! search maximum in neighbourhood 
    309321         ikm1 = MAX(jk-1,1) 
    310322         DO jj = 2, jpjm1 
     
    316328         END DO 
    317329      END DO 
    318       ! large positive value (+zbig) inside land 
     330      !                    ! large positive value (+zbig) inside land 
    319331      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) ) 
    320332      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) ) 
    321       ! search minimum in neighbourhood 
    322       DO jk = 1, jpkm1 
     333      ! 
     334      DO jk = 1, jpkm1     ! search minimum in neighbourhood 
    323335         ikm1 = MAX(jk-1,1) 
    324336         DO jj = 2, jpjm1 
     
    330342         END DO 
    331343      END DO 
    332  
    333       ! restore masked values to zero 
     344      !                    ! restore masked values to zero 
    334345      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) 
    335346      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) 
    336  
    337  
    338       ! 2. Positive and negative part of fluxes and beta terms 
    339       ! ------------------------------------------------------ 
    340  
     347      ! 
     348      ! Positive and negative part of fluxes and beta terms 
     349      ! --------------------------------------------------- 
    341350      DO jk = 1, jpkm1 
    342351         z2dtt = p2dt(jk) 
     
    347356               zneg = MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) ) 
    348357               ! up & down beta terms 
    349                zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt 
     358               zbt = e1e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt 
    350359               zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt 
    351360               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt 
     
    353362         END DO 
    354363      END DO 
     364      ! 
    355365      ! monotonic flux in the k direction, i.e. pcc 
    356366      ! ------------------------------------------- 
     
    366376      END DO 
    367377      ! 
    368       CALL wrk_dealloc( jpi, jpj, jpk, zbetup, zbetdo ) 
     378      CALL wrk_dealloc( jpi,jpj,jpk,  zbetup, zbetdo ) 
    369379      ! 
    370380      IF( nn_timing == 1 )  CALL timing_stop('nonosc_z') 
Note: See TracChangeset for help on using the changeset viewer.