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 13463 for NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/TRA/traldf_iso.F90 – NEMO

Ignore:
Timestamp:
2020-09-14T17:40:34+02:00 (4 years ago)
Author:
andmirek
Message:

Ticket #2195:update to trunk 13461

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

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11351_fldread_with_XIOS

    • Property svn:externals
      •  

        old new  
        33^/utils/build/mk@HEAD         mk 
        44^/utils/tools@HEAD            tools 
        5 ^/vendors/AGRIF/dev@HEAD      ext/AGRIF 
         5^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 
         8 
         9# SETTE 
         10^/utils/CI/sette@13382        sette 
  • NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/TRA/traldf_iso.F90

    r10068 r13463  
    4040 
    4141   !! * Substitutions 
    42 #  include "vectopt_loop_substitute.h90" 
     42#  include "do_loop_substitute.h90" 
     43#  include "domzgr_substitute.h90" 
    4344   !!---------------------------------------------------------------------- 
    4445   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    4849CONTAINS 
    4950 
    50   SUBROUTINE tra_ldf_iso( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   & 
    51       &                                                   pgui, pgvi,   & 
    52       &                                       ptb , ptbb, pta , kjpt, kpass ) 
     51  SUBROUTINE tra_ldf_iso( kt, Kmm, kit000, cdtype, pahu, pahv,                    & 
     52      &                                            pgu , pgv    ,   pgui, pgvi,   & 
     53      &                                       pt , pt2 , pt_rhs , kjpt  , kpass ) 
    5354      !!---------------------------------------------------------------------- 
    5455      !!                  ***  ROUTINE tra_ldf_iso  *** 
     
    8788      !!         difft = 1/(e1e2t*e3t) dk[ zftw ] 
    8889      !!      Add this trend to the general trend (ta,sa): 
    89       !!         pta = pta + difft 
    90       !! 
    91       !! ** Action :   Update pta arrays with the before rotated diffusion 
     90      !!         pt_rhs = pt_rhs + difft 
     91      !! 
     92      !! ** Action :   Update pt_rhs arrays with the before rotated diffusion 
    9293      !!---------------------------------------------------------------------- 
    9394      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index 
     
    9697      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers 
    9798      INTEGER                              , INTENT(in   ) ::   kpass      ! =1/2 first or second passage 
     99      INTEGER                              , INTENT(in   ) ::   Kmm        ! ocean time level index 
    98100      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s] 
    99101      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels 
    100102      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels 
    101       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! tracer (kpass=1) or laplacian of tracer (kpass=2) 
    102       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptbb       ! tracer (only used in kpass=2) 
    103       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend 
     103      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt         ! tracer (kpass=1) or laplacian of tracer (kpass=2) 
     104      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt2        ! tracer (only used in kpass=2) 
     105      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs     ! tracer trend 
    104106      ! 
    105107      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
     
    108110      REAL(wp) ::  zmsku, zahu_w, zabe1, zcof1, zcoef3   ! local scalars 
    109111      REAL(wp) ::  zmskv, zahv_w, zabe2, zcof2, zcoef4   !   -      - 
    110       REAL(wp) ::  zcoef0, ze3w_2, zsign, z2dt, z1_2dt   !   -      - 
     112      REAL(wp) ::  zcoef0, ze3w_2, zsign                 !   -      - 
    111113      REAL(wp), DIMENSION(jpi,jpj)     ::   zdkt, zdk1t, z2d 
    112114      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdit, zdjt, zftu, zftv, ztfw  
     
    124126      l_hst = .FALSE. 
    125127      l_ptr = .FALSE. 
    126       IF( cdtype == 'TRA' .AND. ln_diaptr )                                                 l_ptr = .TRUE.  
     128      IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) )     l_ptr = .TRUE.  
    127129      IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
    128130         &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE. 
    129131      ! 
    130       !                                            ! set time step size (Euler/Leapfrog) 
    131       IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   z2dt =     rdt      ! at nit000   (Euler) 
    132       ELSE                                        ;   z2dt = 2.* rdt      !             (Leapfrog) 
    133       ENDIF 
    134       z1_2dt = 1._wp / z2dt 
    135132      ! 
    136133      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0) 
     
    144141      IF( kpass == 1 ) THEN                  !==  first pass only  ==! 
    145142         ! 
    146          DO jk = 2, jpkm1 
    147             DO jj = 2, jpjm1 
    148                DO ji = fs_2, fs_jpim1   ! vector opt. 
    149                   ! 
    150                   zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          & 
    151                      &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  ) 
    152                   zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          & 
    153                      &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  ) 
    154                      ! 
    155                   zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    & 
    156                      &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku 
    157                   zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    & 
    158                      &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv 
    159                      ! 
    160                   ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   & 
    161                      &               + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) 
    162                END DO 
    163             END DO 
    164          END DO 
     143         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     144            ! 
     145            zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          & 
     146               &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  ) 
     147            zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          & 
     148               &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  ) 
     149               ! 
     150            zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    & 
     151               &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku 
     152            zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    & 
     153               &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv 
     154               ! 
     155            ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   & 
     156               &               + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) 
     157         END_3D 
    165158         ! 
    166159         IF( ln_traldf_msc ) THEN                ! stabilizing vertical diffusivity coefficient 
    167             DO jk = 2, jpkm1 
    168                DO jj = 2, jpjm1 
    169                   DO ji = fs_2, fs_jpim1 
    170                      akz(ji,jj,jk) = 0.25_wp * (                                                                     & 
    171                         &              ( pahu(ji  ,jj,jk) + pahu(ji  ,jj,jk-1) ) / ( e1u(ji  ,jj) * e1u(ji  ,jj) )   & 
    172                         &            + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) )   & 
    173                         &            + ( pahv(ji,jj  ,jk) + pahv(ji,jj  ,jk-1) ) / ( e2v(ji,jj  ) * e2v(ji,jj  ) )   & 
    174                         &            + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) )   ) 
    175                   END DO 
    176                END DO 
    177             END DO 
     160            DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     161               akz(ji,jj,jk) = 0.25_wp * (                                                                     & 
     162                  &              ( pahu(ji  ,jj,jk) + pahu(ji  ,jj,jk-1) ) / ( e1u(ji  ,jj) * e1u(ji  ,jj) )   & 
     163                  &            + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) )   & 
     164                  &            + ( pahv(ji,jj  ,jk) + pahv(ji,jj  ,jk-1) ) / ( e2v(ji,jj  ) * e2v(ji,jj  ) )   & 
     165                  &            + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) )   ) 
     166            END_3D 
    178167            ! 
    179168            IF( ln_traldf_blp ) THEN                ! bilaplacian operator 
    180                DO jk = 2, jpkm1 
    181                   DO jj = 1, jpjm1 
    182                      DO ji = 1, fs_jpim1 
    183                         akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   & 
    184                            &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) )  ) 
    185                      END DO 
    186                   END DO 
    187                END DO 
     169               DO_3D( 1, 0, 1, 0, 2, jpkm1 ) 
     170                  akz(ji,jj,jk) = 16._wp   & 
     171                     &   * ah_wslp2   (ji,jj,jk)   & 
     172                     &   * (  akz     (ji,jj,jk)   & 
     173                     &      + ah_wslp2(ji,jj,jk)   & 
     174                     &        / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) )  ) 
     175               END_3D 
    188176            ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator 
    189                DO jk = 2, jpkm1 
    190                   DO jj = 1, jpjm1 
    191                      DO ji = 1, fs_jpim1 
    192                         ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) 
    193                         zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
    194                         akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt 
    195                      END DO 
    196                   END DO 
    197                END DO 
     177               DO_3D( 1, 0, 1, 0, 2, jpkm1 ) 
     178                  ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 
     179                  zcoef0 = rDt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
     180                  akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * r1_Dt 
     181               END_3D 
    198182           ENDIF 
    199183           ! 
     
    216200 
    217201         ! Horizontal tracer gradient  
    218          DO jk = 1, jpkm1 
    219             DO jj = 1, jpjm1 
    220                DO ji = 1, fs_jpim1   ! vector opt. 
    221                   zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
    222                   zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
    223                END DO 
    224             END DO 
    225          END DO 
     202         DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 
     203            zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
     204            zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
     205         END_3D 
    226206         IF( ln_zps ) THEN      ! botton and surface ocean correction of the horizontal gradient 
    227             DO jj = 1, jpjm1              ! bottom correction (partial bottom cell) 
    228                DO ji = 1, fs_jpim1   ! vector opt. 
    229                   zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)           
    230                   zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 
    231                END DO 
    232             END DO 
     207            DO_2D( 1, 0, 1, 0 ) 
     208               zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)           
     209               zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 
     210            END_2D 
    233211            IF( ln_isfcav ) THEN      ! first wet level beneath a cavity 
    234                DO jj = 1, jpjm1 
    235                   DO ji = 1, fs_jpim1   ! vector opt. 
    236                      IF( miku(ji,jj) > 1 )   zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)           
    237                      IF( mikv(ji,jj) > 1 )   zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)      
    238                   END DO 
    239                END DO 
     212               DO_2D( 1, 0, 1, 0 ) 
     213                  IF( miku(ji,jj) > 1 )   zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)           
     214                  IF( mikv(ji,jj) > 1 )   zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)      
     215               END_2D 
    240216            ENDIF 
    241217         ENDIF 
     
    248224            ! 
    249225            !                             !== Vertical tracer gradient 
    250             zdk1t(:,:) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * wmask(:,:,jk+1)     ! level jk+1 
     226            zdk1t(:,:) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * wmask(:,:,jk+1)     ! level jk+1 
    251227            ! 
    252228            IF( jk == 1 ) THEN   ;   zdkt(:,:) = zdk1t(:,:)                          ! surface: zdkt(jk=1)=zdkt(jk=2) 
    253             ELSE                 ;   zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk) 
     229            ELSE                 ;   zdkt(:,:) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * wmask(:,:,jk) 
    254230            ENDIF 
    255             DO jj = 1 , jpjm1            !==  Horizontal fluxes 
    256                DO ji = 1, fs_jpim1   ! vector opt. 
    257                   zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) 
    258                   zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) 
    259                   ! 
    260                   zmsku = 1. / MAX(  wmask(ji+1,jj,jk  ) + wmask(ji,jj,jk+1)   & 
    261                      &             + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk  ), 1. ) 
    262                   ! 
    263                   zmskv = 1. / MAX(  wmask(ji,jj+1,jk  ) + wmask(ji,jj,jk+1)   & 
    264                      &             + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk  ), 1. ) 
    265                   ! 
    266                   zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku 
    267                   zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv 
    268                   ! 
    269                   zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)   & 
    270                      &               + zcof1 * (  zdkt (ji+1,jj) + zdk1t(ji,jj)      & 
    271                      &                          + zdk1t(ji+1,jj) + zdkt (ji,jj)  )  ) * umask(ji,jj,jk) 
    272                   zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   & 
    273                      &               + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj)      & 
    274                      &                          + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk)                   
    275                END DO 
    276             END DO 
    277             ! 
    278             DO jj = 2 , jpjm1          !== horizontal divergence and add to pta 
    279                DO ji = fs_2, fs_jpim1   ! vector opt. 
    280                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      & 
    281                      &                                           + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   & 
    282                      &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    283                END DO 
    284             END DO 
     231            DO_2D( 1, 0, 1, 0 ) 
     232               zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) 
     233               zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm) 
     234               ! 
     235               zmsku = 1. / MAX(  wmask(ji+1,jj,jk  ) + wmask(ji,jj,jk+1)   & 
     236                  &             + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk  ), 1. ) 
     237               ! 
     238               zmskv = 1. / MAX(  wmask(ji,jj+1,jk  ) + wmask(ji,jj,jk+1)   & 
     239                  &             + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk  ), 1. ) 
     240               ! 
     241               zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku 
     242               zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv 
     243               ! 
     244               zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)   & 
     245                  &               + zcof1 * (  zdkt (ji+1,jj) + zdk1t(ji,jj)      & 
     246                  &                          + zdk1t(ji+1,jj) + zdkt (ji,jj)  )  ) * umask(ji,jj,jk) 
     247               zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   & 
     248                  &               + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj)      & 
     249                  &                          + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk)                   
     250            END_2D 
     251            ! 
     252            DO_2D( 0, 0, 0, 0 ) 
     253               pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    & 
     254                  &       + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   & 
     255                  &                                              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     256            END_2D 
    285257         END DO                                        !   End of slab   
    286258 
     
    288260         !!   III - vertical trend (full) 
    289261         !!---------------------------------------------------------------------- 
    290          ! 
    291          ztfw(1,:,:) = 0._wp     ;     ztfw(jpi,:,:) = 0._wp 
    292262         ! 
    293263         ! Vertical fluxes 
     
    296266         ztfw(:,:, 1 ) = 0._wp      ;      ztfw(:,:,jpk) = 0._wp 
    297267          
    298          DO jk = 2, jpkm1           ! interior (2=<jk=<jpk-1) 
    299             DO jj = 2, jpjm1 
    300                DO ji = fs_2, fs_jpim1   ! vector opt. 
    301                   ! 
    302                   zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          & 
    303                      &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  ) 
    304                   zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          & 
    305                      &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  ) 
    306                      ! 
    307                   zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    & 
    308                      &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku 
    309                   zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    & 
    310                      &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv 
    311                      ! 
    312                   zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk)   !wslpi & j are already w-masked 
    313                   zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) 
    314                   ! 
    315                   ztfw(ji,jj,jk) = zcoef3 * (   zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)      & 
    316                      &                        + zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk)  )   & 
    317                      &           + zcoef4 * (   zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk)      & 
    318                      &                        + zdjt(ji  ,jj-1,jk-1) + zdjt(ji  ,jj  ,jk)  ) 
    319                END DO 
    320             END DO 
    321          END DO 
     268         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     269            ! 
     270            zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          & 
     271               &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  ) 
     272            zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          & 
     273               &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  ) 
     274               ! 
     275            zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    & 
     276               &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku 
     277            zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    & 
     278               &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv 
     279               ! 
     280            zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk)   !wslpi & j are already w-masked 
     281            zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) 
     282            ! 
     283            ztfw(ji,jj,jk) = zcoef3 * (   zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)      & 
     284               &                        + zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk)  )   & 
     285               &           + zcoef4 * (   zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk)      & 
     286               &                        + zdjt(ji  ,jj-1,jk-1) + zdjt(ji  ,jj  ,jk)  ) 
     287         END_3D 
    322288         !                                !==  add the vertical 33 flux  ==! 
    323289         IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz 
    324             DO jk = 2, jpkm1        
    325                DO jj = 1, jpjm1 
    326                   DO ji = fs_2, fs_jpim1 
    327                      ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)   & 
    328                         &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             & 
    329                         &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 
    330                   END DO 
    331                END DO 
    332             END DO 
     290            DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     291               ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)   & 
     292                  &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )               & 
     293                  &                            * (  pt(ji,jj,jk-1,jn) -  pt(ji,jj,jk,jn) ) 
     294            END_3D 
    333295            ! 
    334296         ELSE                                   ! bilaplacian  
    335297            SELECT CASE( kpass ) 
    336298            CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2 
    337                DO jk = 2, jpkm1  
    338                   DO jj = 1, jpjm1 
    339                      DO ji = fs_2, fs_jpim1 
    340                         ztfw(ji,jj,jk) = ztfw(ji,jj,jk)    & 
    341                            &           + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   & 
    342                            &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) 
    343                      END DO 
    344                   END DO 
    345                END DO  
    346             CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb  and ptbb gradients, resp. 
    347                DO jk = 2, jpkm1  
    348                   DO jj = 1, jpjm1 
    349                      DO ji = fs_2, fs_jpim1 
    350                         ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)                      & 
    351                            &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   & 
    352                            &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   ) 
    353                      END DO 
    354                   END DO 
    355                END DO 
     299               DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     300                  ztfw(ji,jj,jk) =   & 
     301                     &  ztfw(ji,jj,jk) + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   & 
     302                     &           * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) 
     303               END_3D 
     304            CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt2 gradients, resp. 
     305               DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     306                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)                  & 
     307                     &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   & 
     308                     &                            +         akz(ji,jj,jk) * ( pt2(ji,jj,jk-1,jn) - pt2(ji,jj,jk,jn) )   ) 
     309               END_3D 
    356310            END SELECT 
    357311         ENDIF 
    358312         !          
    359          DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pta  ==! 
    360             DO jj = 2, jpjm1 
    361                DO ji = fs_2, fs_jpim1   ! vector opt. 
    362                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   & 
    363                      &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    364                END DO 
    365             END DO 
    366          END DO 
     313         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     314            pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  ) * r1_e1e2t(ji,jj)   & 
     315               &                                             / e3t(ji,jj,jk,Kmm) 
     316         END_3D 
    367317         ! 
    368318         IF( ( kpass == 1 .AND. ln_traldf_lap ) .OR.  &     !==  first pass only (  laplacian)  ==! 
Note: See TracChangeset for help on using the changeset viewer.