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 12928 for NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/TRA/traadv_ubs.F90 – NEMO

Ignore:
Timestamp:
2020-05-14T21:46:00+02:00 (4 years ago)
Author:
smueller
Message:

Synchronizing with /NEMO/trunk@12925 (ticket #2170)

Location:
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser

    • Property svn:externals
      •  

        old new  
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 
         8 
         9# SETTE 
         10^/utils/CI/sette@HEAD         sette 
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/TRA/traadv_ubs.F90

    r10425 r12928  
    3838 
    3939   !! * Substitutions 
    40 #  include "vectopt_loop_substitute.h90" 
     40#  include "do_loop_substitute.h90" 
    4141   !!---------------------------------------------------------------------- 
    4242   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    4646CONTAINS 
    4747 
    48    SUBROUTINE tra_adv_ubs( kt, kit000, cdtype, p2dt, pun, pvn, pwn,          & 
    49       &                                                ptb, ptn, pta, kjpt, kn_ubs_v ) 
     48   SUBROUTINE tra_adv_ubs( kt, kit000, cdtype, p2dt, pU, pV, pW,          & 
     49      &                    Kbb, Kmm, pt, kjpt, Krhs, kn_ubs_v ) 
    5050      !!---------------------------------------------------------------------- 
    5151      !!                  ***  ROUTINE tra_adv_ubs  *** 
     
    7777      !!      scheme (kn_ubs_v=4). 
    7878      !! 
    79       !! ** Action : - update pta  with the now advective tracer trends 
     79      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends 
    8080      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T) 
    81       !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) 
     81      !!             - poleward advective heat and salt transport (ln_diaptr=T) 
    8282      !! 
    8383      !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404.  
    8484      !!             Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731Ð1741.  
    8585      !!---------------------------------------------------------------------- 
    86       INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index 
    87       INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index 
    88       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator) 
    89       INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers 
    90       INTEGER                              , INTENT(in   ) ::   kn_ubs_v        ! number of tracers 
    91       REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step 
    92       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean transport components 
    93       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields 
    94       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend  
     86      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index 
     87      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices 
     88      INTEGER                                  , INTENT(in   ) ::   kit000          ! first time step index 
     89      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator) 
     90      INTEGER                                  , INTENT(in   ) ::   kjpt            ! number of tracers 
     91      INTEGER                                  , INTENT(in   ) ::   kn_ubs_v        ! number of tracers 
     92      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step 
     93      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume transport components 
     94      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation 
    9595      ! 
    9696      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
     
    111111      l_ptr = .FALSE. 
    112112      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE. 
    113       IF(   cdtype == 'TRA' .AND. ln_diaptr )                                               l_ptr = .TRUE.  
     113      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )   l_ptr = .TRUE.  
    114114      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
    115115         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE. 
     
    124124         !                                               
    125125         DO jk = 1, jpkm1        !==  horizontal laplacian of before tracer ==! 
    126             DO jj = 1, jpjm1              ! First derivative (masked gradient) 
    127                DO ji = 1, fs_jpim1   ! vector opt. 
    128                   zeeu = e2_e1u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) 
    129                   zeev = e1_e2v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) 
    130                   ztu(ji,jj,jk) = zeeu * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) 
    131                   ztv(ji,jj,jk) = zeev * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) 
    132                END DO 
    133             END DO 
    134             DO jj = 2, jpjm1              ! Second derivative (divergence) 
    135                DO ji = fs_2, fs_jpim1   ! vector opt. 
    136                   zcoef = 1._wp / ( 6._wp * e3t_n(ji,jj,jk) ) 
    137                   zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef 
    138                   zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef 
    139                END DO 
    140             END DO 
     126            DO_2D_10_10 
     127               zeeu = e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) * umask(ji,jj,jk) 
     128               zeev = e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm) * vmask(ji,jj,jk) 
     129               ztu(ji,jj,jk) = zeeu * ( pt(ji+1,jj  ,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) ) 
     130               ztv(ji,jj,jk) = zeev * ( pt(ji  ,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) ) 
     131            END_2D 
     132            DO_2D_00_00 
     133               zcoef = 1._wp / ( 6._wp * e3t(ji,jj,jk,Kmm) ) 
     134               zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef 
     135               zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef 
     136            END_2D 
    141137            !                                     
    142138         END DO          
    143139         CALL lbc_lnk( 'traadv_ubs', zltu, 'T', 1. )   ;    CALL lbc_lnk( 'traadv_ubs', zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn) 
    144140         !     
    145          DO jk = 1, jpkm1        !==  Horizontal advective fluxes  ==!     (UBS) 
    146             DO jj = 1, jpjm1 
    147                DO ji = 1, fs_jpim1   ! vector opt. 
    148                   zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )      ! upstream transport (x2) 
    149                   zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) 
    150                   zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) 
    151                   zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) 
    152                   !                                                  ! 2nd order centered advective fluxes (x2) 
    153                   zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn) ) 
    154                   zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) ) 
    155                   !                                                  ! UBS advective fluxes 
    156                   ztu(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) ) 
    157                   ztv(ji,jj,jk) = 0.5 * ( zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) ) 
    158                END DO 
    159             END DO 
    160          END DO          
    161          ! 
    162          zltu(:,:,:) = pta(:,:,:,jn)      ! store the initial trends before its update 
     141         DO_3D_10_10( 1, jpkm1 ) 
     142            zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) )      ! upstream transport (x2) 
     143            zfm_ui = pU(ji,jj,jk) - ABS( pU(ji,jj,jk) ) 
     144            zfp_vj = pV(ji,jj,jk) + ABS( pV(ji,jj,jk) ) 
     145            zfm_vj = pV(ji,jj,jk) - ABS( pV(ji,jj,jk) ) 
     146            !                                                  ! 2nd order centered advective fluxes (x2) 
     147            zcenut = pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm) ) 
     148            zcenvt = pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm) ) 
     149            !                                                  ! UBS advective fluxes 
     150            ztu(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) ) 
     151            ztv(ji,jj,jk) = 0.5 * ( zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) ) 
     152         END_3D 
     153         ! 
     154         zltu(:,:,:) = pt(:,:,:,jn,Krhs)      ! store the initial trends before its update 
    163155         ! 
    164156         DO jk = 1, jpkm1        !==  add the horizontal advective trend  ==! 
    165             DO jj = 2, jpjm1 
    166                DO ji = fs_2, fs_jpim1   ! vector opt. 
    167                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                        & 
    168                      &             - (  ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk)    & 
    169                      &                + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk)  ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    170                END DO 
    171             END DO 
     157            DO_2D_00_00 
     158               pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs)                        & 
     159                  &             - (  ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk)    & 
     160                  &                + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk)  ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     161            END_2D 
    172162            !                                              
    173163         END DO 
    174164         ! 
    175          zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:)    ! Horizontal advective trend used in vertical 2nd order FCT case 
     165         zltu(:,:,:) = pt(:,:,:,jn,Krhs) - zltu(:,:,:)    ! Horizontal advective trend used in vertical 2nd order FCT case 
    176166         !                                            ! and/or in trend diagnostic (l_trd=T)  
    177167         !                 
    178168         IF( l_trd ) THEN                  ! trend diagnostics 
    179              CALL trd_tra( kt, cdtype, jn, jptra_xad, ztu, pun, ptn(:,:,:,jn) ) 
    180              CALL trd_tra( kt, cdtype, jn, jptra_yad, ztv, pvn, ptn(:,:,:,jn) ) 
     169             CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, ztu, pU, pt(:,:,:,jn,Kmm) ) 
     170             CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, ztv, pV, pt(:,:,:,jn,Kmm) ) 
    181171         END IF 
    182172         !      
     
    193183         CASE(  2  )                   ! 2nd order FCT  
    194184            !          
    195             IF( l_trd )   zltv(:,:,:) = pta(:,:,:,jn)          ! store pta if trend diag. 
     185            IF( l_trd )   zltv(:,:,:) = pt(:,:,:,jn,Krhs)          ! store pt(:,:,:,:,Krhs) if trend diag. 
    196186            ! 
    197187            !                          !*  upstream advection with initial mass fluxes & intermediate update  ==! 
    198             DO jk = 2, jpkm1                 ! Interior value (w-masked) 
    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_wp * (  zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn)  ) * wmask(ji,jj,jk) 
    204                   END DO 
    205                END DO 
    206             END DO  
     188            DO_3D_11_11( 2, jpkm1 ) 
     189               zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) ) 
     190               zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) ) 
     191               ztw(ji,jj,jk) = 0.5_wp * (  zfp_wk * pt(ji,jj,jk,jn,Kbb) + zfm_wk * pt(ji,jj,jk-1,jn,Kbb)  ) * wmask(ji,jj,jk) 
     192            END_3D 
    207193            IF( ln_linssh ) THEN             ! top ocean value (only in linear free surface as ztw has been w-masked) 
    208194               IF( ln_isfcav ) THEN                ! top of the ice-shelf cavities and at the ocean surface 
    209                   DO jj = 1, jpj 
    210                      DO ji = 1, jpi 
    211                         ztw(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
    212                      END DO 
    213                   END DO    
     195                  DO_2D_11_11 
     196                     ztw(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface  
     197                  END_2D 
    214198               ELSE                                ! no cavities: only at the ocean surface 
    215                   ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 
     199                  ztw(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb) 
    216200               ENDIF 
    217201            ENDIF 
    218202            ! 
    219             DO jk = 1, jpkm1           !* trend and after field with monotonic scheme 
    220                DO jj = 2, jpjm1 
    221                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    222                      ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    223                      pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn) +  ztak  
    224                      zti(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + p2dt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk) 
    225                   END DO 
    226                END DO 
    227             END DO 
     203            DO_3D_00_00( 1, jpkm1 ) 
     204               ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     205               pt(ji,jj,jk,jn,Krhs) =   pt(ji,jj,jk,jn,Krhs) +  ztak  
     206               zti(ji,jj,jk)    = ( pt(ji,jj,jk,jn,Kbb) + p2dt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk) 
     207            END_3D 
    228208            CALL lbc_lnk( 'traadv_ubs', zti, 'T', 1. )      ! Lateral boundary conditions on zti, zsi   (unchanged sign) 
    229209            ! 
    230210            !                          !*  anti-diffusive flux : high order minus low order 
    231             DO jk = 2, jpkm1        ! Interior value  (w-masked) 
    232                DO jj = 1, jpj 
    233                   DO ji = 1, jpi 
    234                      ztw(ji,jj,jk) = (   0.5_wp * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )   & 
    235                         &              - ztw(ji,jj,jk)   ) * wmask(ji,jj,jk) 
    236                   END DO 
    237                END DO 
    238             END DO 
     211            DO_3D_11_11( 2, jpkm1 ) 
     212               ztw(ji,jj,jk) = (   0.5_wp * pW(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) )   & 
     213                  &              - ztw(ji,jj,jk)   ) * wmask(ji,jj,jk) 
     214            END_3D 
    239215            !                                            ! top ocean value: high order == upstream  ==>>  zwz=0 
    240216            IF( ln_linssh )   ztw(:,:, 1 ) = 0._wp       ! only ocean surface as interior zwz values have been w-masked 
    241217            ! 
    242             CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt )      !  monotonicity algorithm 
     218            CALL nonosc_z( Kmm, pt(:,:,:,jn,Kbb), ztw, zti, p2dt )      !  monotonicity algorithm 
    243219            ! 
    244220         CASE(  4  )                               ! 4th order COMPACT 
    245             CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )         ! 4th order compact interpolation of T at w-point 
    246             DO jk = 2, jpkm1 
    247                DO jj = 2, jpjm1 
    248                   DO ji = fs_2, fs_jpim1 
    249                      ztw(ji,jj,jk) = pwn(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk) 
    250                   END DO 
    251                END DO 
    252             END DO 
    253             IF( ln_linssh )   ztw(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn)     !!gm ISF & 4th COMPACT doesn't work 
     221            CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw )         ! 4th order compact interpolation of T at w-point 
     222            DO_3D_00_00( 2, jpkm1 ) 
     223               ztw(ji,jj,jk) = pW(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk) 
     224            END_3D 
     225            IF( ln_linssh )   ztw(:,:, 1 ) = pW(:,:,1) * pt(:,:,1,jn,Kmm)     !!gm ISF & 4th COMPACT doesn't work 
    254226            ! 
    255227         END SELECT 
    256228         ! 
    257          DO jk = 1, jpkm1        !  final trend with corrected fluxes 
    258             DO jj = 2, jpjm1  
    259                DO ji = fs_2, fs_jpim1   ! vector opt.    
    260                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    261                END DO 
    262             END DO 
    263          END DO 
     229         DO_3D_00_00( 1, jpkm1 ) 
     230            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     231         END_3D 
    264232         ! 
    265233         IF( l_trd )  THEN       ! vertical advective trend diagnostics 
    266             DO jk = 1, jpkm1                       ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w]) 
    267                DO jj = 2, jpjm1 
    268                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    269                      zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk)                          & 
    270                         &           + ptn(ji,jj,jk,jn) * (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  )   & 
    271                         &                              * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    272                   END DO 
    273                END DO 
    274             END DO 
    275             CALL trd_tra( kt, cdtype, jn, jptra_zad, zltv ) 
     234            DO_3D_00_00( 1, jpkm1 ) 
     235               zltv(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs) - zltv(ji,jj,jk)                          & 
     236                  &           + pt(ji,jj,jk,jn,Kmm) * (  pW(ji,jj,jk) - pW(ji,jj,jk+1)  )   & 
     237                  &                              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     238            END_3D 
     239            CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zltv ) 
    276240         ENDIF 
    277241         ! 
     
    281245 
    282246 
    283    SUBROUTINE nonosc_z( pbef, pcc, paft, p2dt ) 
     247   SUBROUTINE nonosc_z( Kmm, pbef, pcc, paft, p2dt ) 
    284248      !!--------------------------------------------------------------------- 
    285249      !!                    ***  ROUTINE nonosc_z  *** 
     
    294258      !!       in-space based differencing for fluid 
    295259      !!---------------------------------------------------------------------- 
     260      INTEGER , INTENT(in   )                          ::   Kmm    ! time level index 
    296261      REAL(wp), INTENT(in   )                          ::   p2dt   ! tracer time-step 
    297262      REAL(wp),                DIMENSION (jpi,jpj,jpk) ::   pbef   ! before field 
     
    317282      DO jk = 1, jpkm1     ! search maximum in neighbourhood 
    318283         ikm1 = MAX(jk-1,1) 
    319          DO jj = 2, jpjm1 
    320             DO ji = fs_2, fs_jpim1   ! vector opt. 
    321                zbetup(ji,jj,jk) = MAX(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   & 
    322                   &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   & 
    323                   &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  ) 
    324             END DO 
    325          END DO 
     284         DO_2D_00_00 
     285            zbetup(ji,jj,jk) = MAX(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   & 
     286               &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   & 
     287               &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  ) 
     288         END_2D 
    326289      END DO 
    327290      !                    ! large positive value (+zbig) inside land 
     
    331294      DO jk = 1, jpkm1     ! search minimum in neighbourhood 
    332295         ikm1 = MAX(jk-1,1) 
    333          DO jj = 2, jpjm1 
    334             DO ji = fs_2, fs_jpim1   ! vector opt. 
    335                zbetdo(ji,jj,jk) = MIN(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   & 
    336                   &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   & 
    337                   &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  ) 
    338             END DO 
    339          END DO 
     296         DO_2D_00_00 
     297            zbetdo(ji,jj,jk) = MIN(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   & 
     298               &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   & 
     299               &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  ) 
     300         END_2D 
    340301      END DO 
    341302      !                    ! restore masked values to zero 
     
    345306      ! Positive and negative part of fluxes and beta terms 
    346307      ! --------------------------------------------------- 
    347       DO jk = 1, jpkm1 
    348          DO jj = 2, jpjm1 
    349             DO ji = fs_2, fs_jpim1   ! vector opt. 
    350                ! positive & negative part of the flux 
    351                zpos = MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) ) 
    352                zneg = MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) ) 
    353                ! up & down beta terms 
    354                zbt = e1e2t(ji,jj) * e3t_n(ji,jj,jk) / p2dt 
    355                zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt 
    356                zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt 
    357             END DO 
    358          END DO 
    359       END DO 
     308      DO_3D_00_00( 1, jpkm1 ) 
     309         ! positive & negative part of the flux 
     310         zpos = MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) ) 
     311         zneg = MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) ) 
     312         ! up & down beta terms 
     313         zbt = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) / p2dt 
     314         zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt 
     315         zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt 
     316      END_3D 
    360317      ! 
    361318      ! monotonic flux in the k direction, i.e. pcc 
    362319      ! ------------------------------------------- 
    363       DO jk = 2, jpkm1 
    364          DO jj = 2, jpjm1 
    365             DO ji = fs_2, fs_jpim1   ! vector opt. 
    366                za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) ) 
    367                zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) ) 
    368                zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) ) 
    369                pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb ) 
    370             END DO 
    371          END DO 
    372       END DO 
     320      DO_3D_00_00( 2, jpkm1 ) 
     321         za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) ) 
     322         zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) ) 
     323         zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) ) 
     324         pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb ) 
     325      END_3D 
    373326      ! 
    374327   END SUBROUTINE nonosc_z 
Note: See TracChangeset for help on using the changeset viewer.