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 12377 for NEMO/trunk/src/OCE/TRA/traldf_iso.F90 – NEMO

Ignore:
Timestamp:
2020-02-12T15:39:06+01:00 (4 years ago)
Author:
acc
Message:

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

Location:
NEMO/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk

    • 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_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 

  • NEMO/trunk/src/OCE/TRA/traldf_iso.F90

    r11993 r12377  
    4040 
    4141   !! * Substitutions 
    42 #  include "vectopt_loop_substitute.h90" 
     42#  include "do_loop_substitute.h90" 
    4343   !!---------------------------------------------------------------------- 
    4444   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    4848CONTAINS 
    4949 
    50   SUBROUTINE tra_ldf_iso( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   & 
    51       &                                                   pgui, pgvi,   & 
    52       &                                       ptb , ptbb, pta , kjpt, kpass ) 
     50  SUBROUTINE tra_ldf_iso( kt, Kmm, kit000, cdtype, pahu, pahv,                    & 
     51      &                                            pgu , pgv    ,   pgui, pgvi,   & 
     52      &                                       pt , pt2 , pt_rhs , kjpt  , kpass ) 
    5353      !!---------------------------------------------------------------------- 
    5454      !!                  ***  ROUTINE tra_ldf_iso  *** 
    … …  
    8787      !!         difft = 1/(e1e2t*e3t) dk[ zftw ] 
    8888      !!      Add this trend to the general trend (ta,sa): 
    89       !!         pta = pta + difft 
    90       !! 
    91       !! ** Action :   Update pta arrays with the before rotated diffusion 
     89      !!         pt_rhs = pt_rhs + difft 
     90      !! 
     91      !! ** Action :   Update pt_rhs arrays with the before rotated diffusion 
    9292      !!---------------------------------------------------------------------- 
    9393      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index 
    … …  
    9696      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers 
    9797      INTEGER                              , INTENT(in   ) ::   kpass      ! =1/2 first or second passage 
     98      INTEGER                              , INTENT(in   ) ::   Kmm        ! ocean time level index 
    9899      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s] 
    99100      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels 
    100101      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 
     102      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt         ! tracer (kpass=1) or laplacian of tracer (kpass=2) 
     103      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt2        ! tracer (only used in kpass=2) 
     104      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs     ! tracer trend 
    104105      ! 
    105106      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
    … …  
    124125      l_hst = .FALSE. 
    125126      l_ptr = .FALSE. 
    126       IF( cdtype == 'TRA' .AND. ln_diaptr )                                                 l_ptr = .TRUE.  
     127      IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) )     l_ptr = .TRUE.  
    127128      IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
    128129         &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE. 
    … …  
    144145      IF( kpass == 1 ) THEN                  !==  first pass only  ==! 
    145146         ! 
    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 
     147         DO_3D_00_00( 2, jpkm1 ) 
     148            ! 
     149            zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          & 
     150               &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  ) 
     151            zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          & 
     152               &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  ) 
     153               ! 
     154            zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    & 
     155               &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku 
     156            zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    & 
     157               &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv 
     158               ! 
     159            ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   & 
     160               &               + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) 
     161         END_3D 
    165162         ! 
    166163         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 
     164            DO_3D_00_00( 2, jpkm1 ) 
     165               akz(ji,jj,jk) = 0.25_wp * (                                                                     & 
     166                  &              ( pahu(ji  ,jj,jk) + pahu(ji  ,jj,jk-1) ) / ( e1u(ji  ,jj) * e1u(ji  ,jj) )   & 
     167                  &            + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) )   & 
     168                  &            + ( pahv(ji,jj  ,jk) + pahv(ji,jj  ,jk-1) ) / ( e2v(ji,jj  ) * e2v(ji,jj  ) )   & 
     169                  &            + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) )   ) 
     170            END_3D 
    178171            ! 
    179172            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 
     173               DO_3D_10_10( 2, jpkm1 ) 
     174                  akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   & 
     175                     &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) )  ) 
     176               END_3D 
    188177            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 
     178               DO_3D_10_10( 2, jpkm1 ) 
     179                  ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 
     180                  zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
     181                  akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt 
     182               END_3D 
    198183           ENDIF 
    199184           ! 
    … …  
    216201 
    217202         ! 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 
     203         DO_3D_10_10( 1, jpkm1 ) 
     204            zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
     205            zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
     206         END_3D 
    226207         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 
     208            DO_2D_10_10 
     209               zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)           
     210               zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 
     211            END_2D 
    233212            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 
     213               DO_2D_10_10 
     214                  IF( miku(ji,jj) > 1 )   zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)           
     215                  IF( mikv(ji,jj) > 1 )   zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)      
     216               END_2D 
    240217            ENDIF 
    241218         ENDIF 
    … …  
    248225            ! 
    249226            !                             !== Vertical tracer gradient 
    250             zdk1t(:,:) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * wmask(:,:,jk+1)     ! level jk+1 
     227            zdk1t(:,:) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * wmask(:,:,jk+1)     ! level jk+1 
    251228            ! 
    252229            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) 
     230            ELSE                 ;   zdkt(:,:) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * wmask(:,:,jk) 
    254231            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 
     232            DO_2D_10_10 
     233               zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) 
     234               zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm) 
     235               ! 
     236               zmsku = 1. / MAX(  wmask(ji+1,jj,jk  ) + wmask(ji,jj,jk+1)   & 
     237                  &             + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk  ), 1. ) 
     238               ! 
     239               zmskv = 1. / MAX(  wmask(ji,jj+1,jk  ) + wmask(ji,jj,jk+1)   & 
     240                  &             + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk  ), 1. ) 
     241               ! 
     242               zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku 
     243               zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv 
     244               ! 
     245               zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)   & 
     246                  &               + zcof1 * (  zdkt (ji+1,jj) + zdk1t(ji,jj)      & 
     247                  &                          + zdk1t(ji+1,jj) + zdkt (ji,jj)  )  ) * umask(ji,jj,jk) 
     248               zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   & 
     249                  &               + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj)      & 
     250                  &                          + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk)                   
     251            END_2D 
     252            ! 
     253            DO_2D_00_00 
     254               pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      & 
     255                  &                                                 + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   & 
     256                  &                                              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     257            END_2D 
    285258         END DO                                        !   End of slab   
    286259 
    … …  
    288261         !!   III - vertical trend (full) 
    289262         !!---------------------------------------------------------------------- 
    290          ! 
    291          ztfw(fs_2:1,:,:) = 0._wp     ;     ztfw(jpi:fs_jpim1,:,:) = 0._wp   ! avoid to potentially manipulate NaN values 
    292263         ! 
    293264         ! Vertical fluxes 
    … …  
    296267         ztfw(:,:, 1 ) = 0._wp      ;      ztfw(:,:,jpk) = 0._wp 
    297268          
    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 
     269         DO_3D_00_00( 2, jpkm1 ) 
     270            ! 
     271            zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          & 
     272               &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  ) 
     273            zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          & 
     274               &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  ) 
     275               ! 
     276            zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    & 
     277               &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku 
     278            zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    & 
     279               &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv 
     280               ! 
     281            zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk)   !wslpi & j are already w-masked 
     282            zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) 
     283            ! 
     284            ztfw(ji,jj,jk) = zcoef3 * (   zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)      & 
     285               &                        + zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk)  )   & 
     286               &           + zcoef4 * (   zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk)      & 
     287               &                        + zdjt(ji  ,jj-1,jk-1) + zdjt(ji  ,jj  ,jk)  ) 
     288         END_3D 
    322289         !                                !==  add the vertical 33 flux  ==! 
    323290         IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz 
    324             DO jk = 2, jpkm1        
    325                DO jj = 2, 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 
     291            DO_3D_00_00( 2, jpkm1 ) 
     292               ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)   & 
     293                  &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )               & 
     294                  &                            * (  pt(ji,jj,jk-1,jn) -  pt(ji,jj,jk,jn) ) 
     295            END_3D 
    333296            ! 
    334297         ELSE                                   ! bilaplacian  
    335298            SELECT CASE( kpass ) 
    336299            CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2 
    337                DO jk = 2, jpkm1  
    338                   DO jj = 2, 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 = 2, 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 
     300               DO_3D_00_00( 2, jpkm1 ) 
     301                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk)                       & 
     302                     &           + ah_wslp2(ji,jj,jk)  * e1e2t(ji,jj)   & 
     303                     &           * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) 
     304               END_3D 
     305            CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt2 gradients, resp. 
     306               DO_3D_00_00( 2, jpkm1 ) 
     307                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)                  & 
     308                     &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   & 
     309                     &                            +         akz(ji,jj,jk) * ( pt2(ji,jj,jk-1,jn) - pt2(ji,jj,jk,jn) )   ) 
     310               END_3D 
    356311            END SELECT 
    357312         ENDIF 
    358313         !          
    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 
     314         DO_3D_00_00( 1, jpkm1 ) 
     315            pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   & 
     316               &                                              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     317         END_3D 
    367318         ! 
    368319         IF( ( kpass == 1 .AND. ln_traldf_lap ) .OR.  &     !==  first pass only (  laplacian)  ==! 
Note: See TracChangeset for help on using the changeset viewer.