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/DYN/dynzdf.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/DYN/dynzdf.F90

    r11281 r13463  
    3737 
    3838   !! * Substitutions 
    39 #  include "vectopt_loop_substitute.h90" 
     39#  include "do_loop_substitute.h90" 
     40#  include "domzgr_substitute.h90" 
    4041   !!---------------------------------------------------------------------- 
    4142   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    4546CONTAINS 
    4647    
    47    SUBROUTINE dyn_zdf( kt ) 
     48   SUBROUTINE dyn_zdf( kt, Kbb, Kmm, Krhs, puu, pvv, Kaa ) 
    4849      !!---------------------------------------------------------------------- 
    4950      !!                  ***  ROUTINE dyn_zdf  *** 
     
    5455      !! 
    5556      !! ** Method  :  - Leap-Frog time stepping on all trends but the vertical mixing 
    56       !!         ua =         ub + 2*dt *       ua             vector form or linear free surf. 
    57       !!         ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a   otherwise 
     57      !!         u(after) =         u(before) + 2*dt *       u(rhs)                vector form or linear free surf. 
     58      !!         u(after) = ( e3u_b*u(before) + 2*dt * e3u_n*u(rhs) ) / e3u_after   otherwise 
    5859      !!               - update the after velocity with the implicit vertical mixing. 
    5960      !!      This requires to solver the following system:  
    60       !!         ua = ua + 1/e3u_a dk+1[ mi(avm) / e3uw_a dk[ua] ] 
     61      !!         u(after) = u(after) + 1/e3u_after  dk+1[ mi(avm) / e3uw_after dk[ua] ] 
    6162      !!      with the following surface/top/bottom boundary condition: 
    6263      !!      surface: wind stress input (averaged over kt-1/2 & kt+1/2) 
    6364      !!      top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90) 
    6465      !! 
    65       !! ** Action :   (ua,va)   after velocity  
     66      !! ** Action :   (puu(:,:,:,Kaa),pvv(:,:,:,Kaa))   after velocity  
    6667      !!--------------------------------------------------------------------- 
    67       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     68      INTEGER                             , INTENT( in )  ::  kt                  ! ocean time-step index 
     69      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs, Kaa ! ocean time level indices 
     70      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv            ! ocean velocities and RHS of momentum equation 
    6871      ! 
    6972      INTEGER  ::   ji, jj, jk         ! dummy loop indices 
     
    9093         ENDIF 
    9194      ENDIF 
    92       !                             !* set time step 
    93       IF( neuler == 0 .AND. kt == nit000     ) THEN   ;   r2dt =      rdt   ! = rdt (restart with Euler time stepping) 
    94       ELSEIF(               kt <= nit000 + 1 ) THEN   ;   r2dt = 2. * rdt   ! = 2 rdt (leapfrog) 
    95       ENDIF 
    96       ! 
    9795      !                             !* explicit top/bottom drag case 
    98       IF( .NOT.ln_drgimp )   CALL zdf_drg_exp( kt, ub, vb, ua, va )  ! add top/bottom friction trend to (ua,va) 
     96      IF( .NOT.ln_drgimp )   CALL zdf_drg_exp( kt, Kmm, puu(:,:,:,Kbb), pvv(:,:,:,Kbb), puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add top/bottom friction trend to (puu(Kaa),pvv(Kaa)) 
    9997      ! 
    10098      ! 
    10199      IF( l_trddyn )   THEN         !* temporary save of ta and sa trends 
    102100         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) )  
    103          ztrdu(:,:,:) = ua(:,:,:) 
    104          ztrdv(:,:,:) = va(:,:,:) 
    105       ENDIF 
    106       ! 
    107       !              !==  RHS: Leap-Frog time stepping on all trends but the vertical mixing  ==!   (put in ua,va) 
     101         ztrdu(:,:,:) = puu(:,:,:,Krhs) 
     102         ztrdv(:,:,:) = pvv(:,:,:,Krhs) 
     103      ENDIF 
     104      ! 
     105      !              !==  RHS: Leap-Frog time stepping on all trends but the vertical mixing  ==!   (put in puu(:,:,:,Kaa),pvv(:,:,:,Kaa)) 
    108106      ! 
    109107      !                    ! time stepping except vertical diffusion 
    110108      IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity 
    111          DO jk = 1, jpkm1 
    112             ua(:,:,jk) = ( ub(:,:,jk) + r2dt * ua(:,:,jk) ) * umask(:,:,jk) 
    113             va(:,:,jk) = ( vb(:,:,jk) + r2dt * va(:,:,jk) ) * vmask(:,:,jk) 
    114          END DO 
     109         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     110            puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kbb) + rDt * puu(ji,jj,jk,Krhs) ) * umask(ji,jj,jk) 
     111            pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kbb) + rDt * pvv(ji,jj,jk,Krhs) ) * vmask(ji,jj,jk) 
     112         END_3D 
    115113      ELSE                                      ! applied on thickness weighted velocity 
    116          DO jk = 1, jpkm1 
    117             ua(:,:,jk) = (         e3u_b(:,:,jk) * ub(:,:,jk)  & 
    118                &          + r2dt * e3u_n(:,:,jk) * ua(:,:,jk)  ) / e3u_a(:,:,jk) * umask(:,:,jk) 
    119             va(:,:,jk) = (         e3v_b(:,:,jk) * vb(:,:,jk)  & 
    120                &          + r2dt * e3v_n(:,:,jk) * va(:,:,jk)  ) / e3v_a(:,:,jk) * vmask(:,:,jk) 
    121          END DO 
     114         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     115            puu(ji,jj,jk,Kaa) = (         e3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb )  & 
     116               &                  + rDt * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Krhs)  ) & 
     117               &                        / e3u(ji,jj,jk,Kaa) * umask(ji,jj,jk) 
     118            pvv(ji,jj,jk,Kaa) = (         e3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb )  & 
     119               &                  + rDt * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Krhs)  ) & 
     120               &                        / e3v(ji,jj,jk,Kaa) * vmask(ji,jj,jk) 
     121         END_3D 
    122122      ENDIF 
    123123      !                    ! add top/bottom friction  
     
    125125      !     J. Chanut: The bottom stress is computed considering after barotropic velocities, which does  
    126126      !                not lead to the effective stress seen over the whole barotropic loop.  
    127       !     G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a 
     127      !     G. Madec : in linear free surface, e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) = e3u_0, so systematic use of e3u(:,:,:,Kaa) 
    128128      IF( ln_drgimp .AND. ln_dynspg_ts ) THEN 
    129          DO jk = 1, jpkm1        ! remove barotropic velocities 
    130             ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk) 
    131             va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk) 
    132          END DO 
    133          DO jj = 2, jpjm1        ! Add bottom/top stress due to barotropic component only 
    134             DO ji = fs_2, fs_jpim1   ! vector opt. 
    135                iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
    136                ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
    137                ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) 
    138                ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) 
    139                ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua 
    140                va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va 
    141             END DO 
    142          END DO 
     129         DO_3D( 0, 0, 0, 0, 1, jpkm1 )      ! remove barotropic velocities 
     130            puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - uu_b(ji,jj,Kaa) ) * umask(ji,jj,jk) 
     131            pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - vv_b(ji,jj,Kaa) ) * vmask(ji,jj,jk) 
     132         END_3D 
     133         DO_2D( 0, 0, 0, 0 ) 
     134            iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
     135            ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
     136            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm)    & 
     137               &             + r_vvl   * e3u(ji,jj,iku,Kaa) 
     138            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    & 
     139               &             + r_vvl   * e3v(ji,jj,ikv,Kaa) 
     140            puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua 
     141            pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va 
     142         END_2D 
    143143         IF( ln_isfcav ) THEN    ! Ocean cavities (ISF) 
    144             DO jj = 2, jpjm1         
    145                DO ji = fs_2, fs_jpim1   ! vector opt. 
    146                   iku = miku(ji,jj)         ! top ocean level at u- and v-points  
    147                   ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points) 
    148                   ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) 
    149                   ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) 
    150                   ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua 
    151                   va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va 
    152                END DO 
    153             END DO 
     144            DO_2D( 0, 0, 0, 0 ) 
     145               iku = miku(ji,jj)         ! top ocean level at u- and v-points  
     146               ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points) 
     147               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm)    & 
     148                  &             + r_vvl   * e3u(ji,jj,iku,Kaa) 
     149               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    & 
     150                  &             + r_vvl   * e3v(ji,jj,ikv,Kaa) 
     151               puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + rDt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua 
     152               pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va 
     153            END_2D 
    154154         END IF 
    155155      ENDIF 
     
    158158      ! 
    159159      !                    !* Matrix construction 
    160       zdt = r2dt * 0.5 
     160      zdt = rDt * 0.5 
    161161      IF( ln_zad_Aimp ) THEN   !! 
    162162         SELECT CASE( nldf_dyn ) 
    163163         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu) 
    164             DO jk = 1, jpkm1 
    165                DO jj = 2, jpjm1  
    166                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    167                      ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point 
    168                      zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
    169                         &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
    170                      zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
    171                         &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
    172                      zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua 
    173                      zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua 
    174                      zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp )  
    175                      zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp ) 
    176                      zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 
    177                   END DO 
    178                END DO 
    179             END DO 
     164            DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     165               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    & 
     166                  &             + r_vvl   * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point 
     167               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
     168                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  ) 
     169               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
     170                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 
     171               zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua 
     172               zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua 
     173               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp )  
     174               zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp ) 
     175               zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 
     176            END_3D 
    180177         CASE DEFAULT               ! iso-level lateral mixing 
    181             DO jk = 1, jpkm1 
    182                DO jj = 2, jpjm1  
    183                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    184                      ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point 
    185                      zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
    186                      zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
    187                      zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua 
    188                      zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua 
    189                      zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 
    190                      zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp ) 
    191                      zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 
    192                   END DO 
    193                END DO 
    194             END DO 
     178            DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     179               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    &    ! after scale factor at U-point 
     180                  &             + r_vvl   * e3u(ji,jj,jk,Kaa) 
     181               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) )   & 
     182                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  ) 
     183               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) )   & 
     184                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 
     185               zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua 
     186               zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua 
     187               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 
     188               zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp ) 
     189               zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 
     190            END_3D 
    195191         END SELECT 
    196          DO jj = 2, jpjm1     !* Surface boundary conditions 
    197             DO ji = fs_2, fs_jpim1   ! vector opt. 
    198                zwi(ji,jj,1) = 0._wp 
    199                ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 
    200                zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) ) / ( ze3ua * e3uw_n(ji,jj,2) ) * wumask(ji,jj,2) 
    201                zWus = ( wi(ji  ,jj,2) +  wi(ji+1,jj,2) ) / ze3ua 
    202                zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp ) 
    203                zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ) ) 
    204             END DO 
    205          END DO 
     192         DO_2D( 0, 0, 0, 0 ) 
     193            zwi(ji,jj,1) = 0._wp 
     194            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm)    & 
     195               &             + r_vvl   * e3u(ji,jj,1,Kaa) 
     196            zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) )   & 
     197               &         / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2) 
     198            zWus = ( wi(ji  ,jj,2) +  wi(ji+1,jj,2) ) / ze3ua 
     199            zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp ) 
     200            zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ) ) 
     201         END_2D 
    206202      ELSE 
    207203         SELECT CASE( nldf_dyn ) 
    208204         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu) 
    209             DO jk = 1, jpkm1 
    210                DO jj = 2, jpjm1  
    211                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    212                      ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point 
    213                      zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
    214                         &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
    215                      zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
    216                         &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
    217                      zwi(ji,jj,jk) = zzwi 
    218                      zws(ji,jj,jk) = zzws 
    219                      zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
    220                   END DO 
    221                END DO 
    222             END DO 
     205            DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     206               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    & 
     207                  &             + r_vvl   * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point 
     208               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
     209                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  ) 
     210               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
     211                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 
     212               zwi(ji,jj,jk) = zzwi 
     213               zws(ji,jj,jk) = zzws 
     214               zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     215            END_3D 
    223216         CASE DEFAULT               ! iso-level lateral mixing 
    224             DO jk = 1, jpkm1 
    225                DO jj = 2, jpjm1  
    226                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    227                      ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point 
    228                      zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
    229                      zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
    230                      zwi(ji,jj,jk) = zzwi 
    231                      zws(ji,jj,jk) = zzws 
    232                      zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
    233                   END DO 
    234                END DO 
    235             END DO 
     217            DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     218               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    & 
     219                  &             + r_vvl   * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point 
     220               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) )    & 
     221                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  ) 
     222               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) )    & 
     223                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 
     224               zwi(ji,jj,jk) = zzwi 
     225               zws(ji,jj,jk) = zzws 
     226               zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     227            END_3D 
    236228         END SELECT 
    237          DO jj = 2, jpjm1     !* Surface boundary conditions 
    238             DO ji = fs_2, fs_jpim1   ! vector opt. 
    239                zwi(ji,jj,1) = 0._wp 
    240                zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
    241             END DO 
    242          END DO 
     229         DO_2D( 0, 0, 0, 0 ) 
     230            zwi(ji,jj,1) = 0._wp 
     231            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
     232         END_2D 
    243233      ENDIF 
    244234      ! 
     
    251241      ! 
    252242      IF ( ln_drgimp ) THEN      ! implicit bottom friction 
    253          DO jj = 2, jpjm1 
    254             DO ji = 2, jpim1 
    255                iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points 
    256                ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point 
    257                zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 
    258             END DO 
    259          END DO 
     243         DO_2D( 0, 0, 0, 0 ) 
     244            iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points 
     245            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm)    & 
     246               &             + r_vvl   * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point 
     247            zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 
     248         END_2D 
    260249         IF ( ln_isfcav ) THEN   ! top friction (always implicit) 
    261             DO jj = 2, jpjm1 
    262                DO ji = 2, jpim1 
    263                   !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed 
    264                   iku = miku(ji,jj)       ! ocean top level at u- and v-points  
    265                   ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point 
    266                   zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 
    267                END DO 
    268             END DO 
     250            DO_2D( 0, 0, 0, 0 ) 
     251               !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed 
     252               iku = miku(ji,jj)       ! ocean top level at u- and v-points  
     253               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm)    & 
     254                  &             + r_vvl   * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point 
     255               zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 
     256            END_2D 
    269257         END IF 
    270258      ENDIF 
     
    282270      !   m is decomposed in the product of an upper and a lower triangular matrix 
    283271      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 
    284       !   The solution (the after velocity) is in ua 
     272      !   The solution (the after velocity) is in puu(:,:,:,Kaa) 
    285273      !----------------------------------------------------------------------- 
    286274      ! 
    287       DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    288          DO jj = 2, jpjm1    
    289             DO ji = fs_2, fs_jpim1   ! vector opt. 
    290                zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
    291             END DO 
    292          END DO 
    293       END DO 
    294       ! 
    295       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==! 
    296          DO ji = fs_2, fs_jpim1   ! vector opt. 
    297             ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1)  
    298             ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
    299                &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1)  
    300          END DO 
    301       END DO 
    302       DO jk = 2, jpkm1 
    303          DO jj = 2, jpjm1 
    304             DO ji = fs_2, fs_jpim1 
    305                ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 
    306             END DO 
    307          END DO 
    308       END DO 
    309       ! 
    310       DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==! 
    311          DO ji = fs_2, fs_jpim1   ! vector opt. 
    312             ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
    313          END DO 
    314       END DO 
    315       DO jk = jpk-2, 1, -1 
    316          DO jj = 2, jpjm1 
    317             DO ji = fs_2, fs_jpim1 
    318                ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
    319             END DO 
    320          END DO 
    321       END DO 
     275      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     276         zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
     277      END_3D 
     278      ! 
     279      DO_2D( 0, 0, 0, 0 ) 
     280         ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm)    & 
     281            &             + r_vvl   * e3u(ji,jj,1,Kaa)  
     282         puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
     283            &                                      / ( ze3ua * rho0 ) * umask(ji,jj,1)  
     284      END_2D 
     285      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     286         puu(ji,jj,jk,Kaa) = puu(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * puu(ji,jj,jk-1,Kaa) 
     287      END_3D 
     288      ! 
     289      DO_2D( 0, 0, 0, 0 ) 
     290         puu(ji,jj,jpkm1,Kaa) = puu(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1) 
     291      END_2D 
     292      DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 ) 
     293         puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - zws(ji,jj,jk) * puu(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk) 
     294      END_3D 
    322295      ! 
    323296      !              !==  Vertical diffusion on v  ==! 
    324297      ! 
    325298      !                       !* Matrix construction 
    326       zdt = r2dt * 0.5 
     299      zdt = rDt * 0.5 
    327300      IF( ln_zad_Aimp ) THEN   !! 
    328301         SELECT CASE( nldf_dyn ) 
    329302         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzv) 
    330             DO jk = 1, jpkm1 
    331                DO jj = 2, jpjm1  
    332                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    333                      ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point 
    334                      zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
    335                         &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
    336                      zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
    337                         &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
    338                      zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va 
    339                      zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va 
    340                      zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp ) 
    341                      zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp ) 
    342                      zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 
    343                   END DO 
    344                END DO 
    345             END DO 
     303            DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     304               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    & 
     305                  &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point 
     306               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
     307                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  ) 
     308               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
     309                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 
     310               zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va 
     311               zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va 
     312               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp ) 
     313               zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp ) 
     314               zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 
     315            END_3D 
    346316         CASE DEFAULT               ! iso-level lateral mixing 
    347             DO jk = 1, jpkm1 
    348                DO jj = 2, jpjm1  
    349                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    350                      ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point 
    351                      zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
    352                      zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
    353                      zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va 
    354                      zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va 
    355                      zwi(ji,jj,jk) = zzwi  + zdt * MIN( zWvi, 0._wp ) 
    356                      zws(ji,jj,jk) = zzws  - zdt * MAX( zWvs, 0._wp ) 
    357                      zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 
    358                   END DO 
    359                END DO 
    360             END DO 
     317            DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     318               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    & 
     319                  &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point 
     320               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) )    & 
     321                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  ) 
     322               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) )    & 
     323                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 
     324               zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va 
     325               zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va 
     326               zwi(ji,jj,jk) = zzwi  + zdt * MIN( zWvi, 0._wp ) 
     327               zws(ji,jj,jk) = zzws  - zdt * MAX( zWvs, 0._wp ) 
     328               zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 
     329            END_3D 
    361330         END SELECT 
    362          DO jj = 2, jpjm1     !* Surface boundary conditions 
    363             DO ji = fs_2, fs_jpim1   ! vector opt. 
    364                zwi(ji,jj,1) = 0._wp 
    365                ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 
    366                zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw_n(ji,jj,2) ) * wvmask(ji,jj,2) 
    367                zWvs = ( wi(ji,jj  ,2) +  wi(ji,jj+1,2) ) / ze3va 
    368                zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp ) 
    369                zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) ) 
    370             END DO 
    371          END DO 
     331         DO_2D( 0, 0, 0, 0 ) 
     332            zwi(ji,jj,1) = 0._wp 
     333            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm)    & 
     334               &             + r_vvl   * e3v(ji,jj,1,Kaa) 
     335            zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) )    & 
     336               &         / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2) 
     337            zWvs = ( wi(ji,jj  ,2) +  wi(ji,jj+1,2) ) / ze3va 
     338            zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp ) 
     339            zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) ) 
     340         END_2D 
    372341      ELSE 
    373342         SELECT CASE( nldf_dyn ) 
    374343         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu) 
    375             DO jk = 1, jpkm1 
    376                DO jj = 2, jpjm1    
    377                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    378                      ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point 
    379                      zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
    380                         &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
    381                      zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
    382                         &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
    383                      zwi(ji,jj,jk) = zzwi 
    384                      zws(ji,jj,jk) = zzws 
    385                      zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
    386                   END DO 
    387                END DO 
    388             END DO 
     344            DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     345               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    & 
     346                  &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point 
     347               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
     348                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  ) 
     349               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
     350                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 
     351               zwi(ji,jj,jk) = zzwi 
     352               zws(ji,jj,jk) = zzws 
     353               zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     354            END_3D 
    389355         CASE DEFAULT               ! iso-level lateral mixing 
    390             DO jk = 1, jpkm1 
    391                DO jj = 2, jpjm1    
    392                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    393                      ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point 
    394                      zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
    395                      zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
    396                      zwi(ji,jj,jk) = zzwi 
    397                      zws(ji,jj,jk) = zzws 
    398                      zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
    399                   END DO 
    400                END DO 
    401             END DO 
     356            DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     357               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    & 
     358                  &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point 
     359               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) )    & 
     360                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  ) 
     361               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) )    & 
     362                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 
     363               zwi(ji,jj,jk) = zzwi 
     364               zws(ji,jj,jk) = zzws 
     365               zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     366            END_3D 
    402367         END SELECT 
    403          DO jj = 2, jpjm1        !* Surface boundary conditions 
    404             DO ji = fs_2, fs_jpim1   ! vector opt. 
    405                zwi(ji,jj,1) = 0._wp 
    406                zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
    407             END DO 
    408          END DO 
     368         DO_2D( 0, 0, 0, 0 ) 
     369            zwi(ji,jj,1) = 0._wp 
     370            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
     371         END_2D 
    409372      ENDIF 
    410373      ! 
     
    416379      ! 
    417380      IF( ln_drgimp ) THEN 
    418          DO jj = 2, jpjm1 
    419             DO ji = 2, jpim1 
    420                ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points) 
    421                ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point 
    422                zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va            
    423             END DO 
    424          END DO 
     381         DO_2D( 0, 0, 0, 0 ) 
     382            ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points) 
     383            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    & 
     384               &             + r_vvl   * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point 
     385            zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va            
     386         END_2D 
    425387         IF ( ln_isfcav ) THEN 
    426             DO jj = 2, jpjm1 
    427                DO ji = 2, jpim1 
    428                   ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points) 
    429                   ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point 
    430                   zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va 
    431                END DO 
    432             END DO 
     388            DO_2D( 0, 0, 0, 0 ) 
     389               ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points) 
     390               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    & 
     391                  &             + r_vvl   * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point 
     392               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va 
     393            END_2D 
    433394         ENDIF 
    434395      ENDIF 
     
    449410      !----------------------------------------------------------------------- 
    450411      ! 
    451       DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    452          DO jj = 2, jpjm1    
    453             DO ji = fs_2, fs_jpim1   ! vector opt. 
    454                zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
    455             END DO 
    456          END DO 
    457       END DO 
    458       ! 
    459       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==! 
    460          DO ji = fs_2, fs_jpim1   ! vector opt.           
    461             ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1)  
    462             va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
    463                &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1)  
    464          END DO 
    465       END DO 
    466       DO jk = 2, jpkm1 
    467          DO jj = 2, jpjm1 
    468             DO ji = fs_2, fs_jpim1   ! vector opt. 
    469                va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 
    470             END DO 
    471          END DO 
    472       END DO 
    473       ! 
    474       DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==! 
    475          DO ji = fs_2, fs_jpim1   ! vector opt. 
    476             va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
    477          END DO 
    478       END DO 
    479       DO jk = jpk-2, 1, -1 
    480          DO jj = 2, jpjm1 
    481             DO ji = fs_2, fs_jpim1 
    482                va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
    483             END DO 
    484          END DO 
    485       END DO 
     412      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     413         zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
     414      END_3D 
     415      ! 
     416      DO_2D( 0, 0, 0, 0 ) 
     417         ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm)    & 
     418            &             + r_vvl   * e3v(ji,jj,1,Kaa)  
     419         pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
     420            &                                      / ( ze3va * rho0 ) * vmask(ji,jj,1)  
     421      END_2D 
     422      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     423         pvv(ji,jj,jk,Kaa) = pvv(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * pvv(ji,jj,jk-1,Kaa) 
     424      END_3D 
     425      ! 
     426      DO_2D( 0, 0, 0, 0 ) 
     427         pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1) 
     428      END_2D 
     429      DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 ) 
     430         pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - zws(ji,jj,jk) * pvv(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk) 
     431      END_3D 
    486432      ! 
    487433      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics 
    488          ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:) 
    489          ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:) 
    490          CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt ) 
     434         ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) / rDt - ztrdu(:,:,:) 
     435         ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) / rDt - ztrdv(:,:,:) 
     436         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt, Kmm ) 
    491437         DEALLOCATE( ztrdu, ztrdv )  
    492438      ENDIF 
    493439      !                                          ! print mean trends (used for debugging) 
    494       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask,               & 
    495          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     440      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' zdf  - Ua: ', mask1=umask,               & 
     441         &                                  tab3d_2=pvv(:,:,:,Kaa), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    496442         ! 
    497443      IF( ln_timing )   CALL timing_stop('dyn_zdf') 
Note: See TracChangeset for help on using the changeset viewer.