New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 12928 for NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DYN/dynzdf.F90 – NEMO

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

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

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

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser

    • Property svn:externals
      •  

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

    r12178 r12928  
    3737 
    3838   !! * Substitutions 
    39 #  include "vectopt_loop_substitute.h90" 
     39#  include "do_loop_substitute.h90" 
    4040   !!---------------------------------------------------------------------- 
    4141   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    4545CONTAINS 
    4646    
    47    SUBROUTINE dyn_zdf( kt ) 
     47   SUBROUTINE dyn_zdf( kt, Kbb, Kmm, Krhs, puu, pvv, Kaa ) 
    4848      !!---------------------------------------------------------------------- 
    4949      !!                  ***  ROUTINE dyn_zdf  *** 
     
    5454      !! 
    5555      !! ** 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 
     56      !!         u(after) =         u(before) + 2*dt *       u(rhs)                vector form or linear free surf. 
     57      !!         u(after) = ( e3u_b*u(before) + 2*dt * e3u_n*u(rhs) ) / e3u(after)   otherwise 
    5858      !!               - update the after velocity with the implicit vertical mixing. 
    5959      !!      This requires to solver the following system:  
    60       !!         ua = ua + 1/e3u_a dk+1[ mi(avm) / e3uw_a dk[ua] ] 
     60      !!         u(after) = u(after) + 1/e3u(after) dk+1[ mi(avm) / e3uw(after) dk[ua] ] 
    6161      !!      with the following surface/top/bottom boundary condition: 
    6262      !!      surface: wind stress input (averaged over kt-1/2 & kt+1/2) 
    6363      !!      top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90) 
    6464      !! 
    65       !! ** Action :   (ua,va)   after velocity  
     65      !! ** Action :   (puu(:,:,:,Kaa),pvv(:,:,:,Kaa))   after velocity  
    6666      !!--------------------------------------------------------------------- 
    67       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     67      INTEGER                             , INTENT( in )  ::  kt                  ! ocean time-step index 
     68      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs, Kaa ! ocean time level indices 
     69      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv            ! ocean velocities and RHS of momentum equation 
    6870      ! 
    6971      INTEGER  ::   ji, jj, jk         ! dummy loop indices 
     
    9092         ENDIF 
    9193      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       ! 
    9794      !                             !* 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) 
     95      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)) 
    9996      ! 
    10097      ! 
    10198      IF( l_trddyn )   THEN         !* temporary save of ta and sa trends 
    10299         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) 
     100         ztrdu(:,:,:) = puu(:,:,:,Krhs) 
     101         ztrdv(:,:,:) = pvv(:,:,:,Krhs) 
     102      ENDIF 
     103      ! 
     104      !              !==  RHS: Leap-Frog time stepping on all trends but the vertical mixing  ==!   (put in puu(:,:,:,Kaa),pvv(:,:,:,Kaa)) 
    108105      ! 
    109106      !                    ! time stepping except vertical diffusion 
    110107      IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity 
    111108         DO jk = 1, jpkm1 
    112             ua(:,:,jk) = ( ub(:,:,jk) + r2dt * ua(:,:,jk) ) * umask(:,:,jk) 
    113             va(:,:,jk) = ( vb(:,:,jk) + r2dt * va(:,:,jk) ) * vmask(:,:,jk) 
     109            puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kbb) + rDt * puu(:,:,jk,Krhs) ) * umask(:,:,jk) 
     110            pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kbb) + rDt * pvv(:,:,jk,Krhs) ) * vmask(:,:,jk) 
    114111         END DO 
    115112      ELSE                                      ! applied on thickness weighted velocity 
    116113         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) 
     114            puu(:,:,jk,Kaa) = (         e3u(:,:,jk,Kbb) * puu(:,:,jk,Kbb)  & 
     115               &          + rDt * e3u(:,:,jk,Kmm) * puu(:,:,jk,Krhs)  ) / e3u(:,:,jk,Kaa) * umask(:,:,jk) 
     116            pvv(:,:,jk,Kaa) = (         e3v(:,:,jk,Kbb) * pvv(:,:,jk,Kbb)  & 
     117               &          + rDt * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Krhs)  ) / e3v(:,:,jk,Kaa) * vmask(:,:,jk) 
    121118         END DO 
    122119      ENDIF 
     
    125122      !     J. Chanut: The bottom stress is computed considering after barotropic velocities, which does  
    126123      !                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 
     124      !     G. Madec : in linear free surface, e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) = e3u_0, so systematic use of e3u(:,:,:,Kaa) 
    128125      IF( ln_drgimp .AND. ln_dynspg_ts ) THEN 
    129126         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) 
     127            puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - uu_b(:,:,Kaa) ) * umask(:,:,jk) 
     128            pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - vv_b(:,:,Kaa) ) * vmask(:,:,jk) 
    132129         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 
     130         DO_2D_00_00 
     131            iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
     132            ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
     133            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) 
     134            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) 
     135            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 
     136            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 
     137         END_2D 
    143138         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 
     139            DO_2D_00_00 
     140               iku = miku(ji,jj)         ! top ocean level at u- and v-points  
     141               ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points) 
     142               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) 
     143               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) 
     144               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 
     145               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 
     146            END_2D 
    154147         END IF 
    155148      ENDIF 
     
    158151      ! 
    159152      !                    !* Matrix construction 
    160       zdt = r2dt * 0.5 
     153      zdt = rDt * 0.5 
    161154      IF( ln_zad_Aimp ) THEN   !! 
    162155         SELECT CASE( nldf_dyn ) 
    163156         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 
     157            DO_3D_00_00( 1, jpkm1 ) 
     158               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point 
     159               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
     160                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  ) 
     161               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
     162                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 
     163               zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua 
     164               zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua 
     165               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp )  
     166               zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp ) 
     167               zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 
     168            END_3D 
    180169         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 
     170            DO_3D_00_00( 1, jpkm1 ) 
     171               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point 
     172               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  ) 
     173               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 
     174               zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua 
     175               zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua 
     176               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 
     177               zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp ) 
     178               zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 
     179            END_3D 
    195180         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 
     181         DO_2D_00_00 
     182            zwi(ji,jj,1) = 0._wp 
     183            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa) 
     184            zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) ) / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2) 
     185            zWus = ( wi(ji  ,jj,2) +  wi(ji+1,jj,2) ) / ze3ua 
     186            zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp ) 
     187            zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ) ) 
     188         END_2D 
    206189      ELSE 
    207190         SELECT CASE( nldf_dyn ) 
    208191         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 
     192            DO_3D_00_00( 1, jpkm1 ) 
     193               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point 
     194               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
     195                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  ) 
     196               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
     197                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 
     198               zwi(ji,jj,jk) = zzwi 
     199               zws(ji,jj,jk) = zzws 
     200               zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     201            END_3D 
    223202         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 
     203            DO_3D_00_00( 1, jpkm1 ) 
     204               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point 
     205               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  ) 
     206               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 
     207               zwi(ji,jj,jk) = zzwi 
     208               zws(ji,jj,jk) = zzws 
     209               zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     210            END_3D 
    236211         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 
     212         DO_2D_00_00 
     213            zwi(ji,jj,1) = 0._wp 
     214            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
     215         END_2D 
    243216      ENDIF 
    244217      ! 
     
    251224      ! 
    252225      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 
     226         DO_2D_00_00 
     227            iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points 
     228            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point 
     229            zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 
     230         END_2D 
    260231         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 
     232            DO_2D_00_00 
     233               !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed 
     234               iku = miku(ji,jj)       ! ocean top level at u- and v-points  
     235               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point 
     236               zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 
     237            END_2D 
    269238         END IF 
    270239      ENDIF 
     
    282251      !   m is decomposed in the product of an upper and a lower triangular matrix 
    283252      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 
    284       !   The solution (the after velocity) is in ua 
     253      !   The solution (the after velocity) is in puu(:,:,:,Kaa) 
    285254      !----------------------------------------------------------------------- 
    286255      ! 
    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 
     256      DO_3D_00_00( 2, jpkm1 ) 
     257         zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
     258      END_3D 
     259      ! 
     260      DO_2D_00_00 
     261         ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa)  
     262         puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
     263            &                                      / ( ze3ua * rho0 ) * umask(ji,jj,1)  
     264      END_2D 
     265      DO_3D_00_00( 2, jpkm1 ) 
     266         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) 
     267      END_3D 
     268      ! 
     269      DO_2D_00_00 
     270         puu(ji,jj,jpkm1,Kaa) = puu(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1) 
     271      END_2D 
     272      DO_3DS_00_00( jpk-2, 1, -1 ) 
     273         puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - zws(ji,jj,jk) * puu(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk) 
     274      END_3D 
    322275      ! 
    323276      !              !==  Vertical diffusion on v  ==! 
    324277      ! 
    325278      !                       !* Matrix construction 
    326       zdt = r2dt * 0.5 
     279      zdt = rDt * 0.5 
    327280      IF( ln_zad_Aimp ) THEN   !! 
    328281         SELECT CASE( nldf_dyn ) 
    329282         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 
     283            DO_3D_00_00( 1, jpkm1 ) 
     284               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point 
     285               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
     286                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  ) 
     287               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
     288                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 
     289               zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va 
     290               zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va 
     291               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp ) 
     292               zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp ) 
     293               zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 
     294            END_3D 
    346295         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 
     296            DO_3D_00_00( 1, jpkm1 ) 
     297               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point 
     298               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  ) 
     299               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 
     300               zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va 
     301               zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va 
     302               zwi(ji,jj,jk) = zzwi  + zdt * MIN( zWvi, 0._wp ) 
     303               zws(ji,jj,jk) = zzws  - zdt * MAX( zWvs, 0._wp ) 
     304               zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 
     305            END_3D 
    361306         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 
     307         DO_2D_00_00 
     308            zwi(ji,jj,1) = 0._wp 
     309            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa) 
     310            zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2) 
     311            zWvs = ( wi(ji,jj  ,2) +  wi(ji,jj+1,2) ) / ze3va 
     312            zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp ) 
     313            zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) ) 
     314         END_2D 
    372315      ELSE 
    373316         SELECT CASE( nldf_dyn ) 
    374317         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 
     318            DO_3D_00_00( 1, jpkm1 ) 
     319               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + 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  ) + akzv(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) + akzv(ji,jj,jk+1) )   & 
     323                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 
     324               zwi(ji,jj,jk) = zzwi 
     325               zws(ji,jj,jk) = zzws 
     326               zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     327            END_3D 
    389328         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 
     329            DO_3D_00_00( 1, jpkm1 ) 
     330               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point 
     331               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  ) 
     332               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 
     333               zwi(ji,jj,jk) = zzwi 
     334               zws(ji,jj,jk) = zzws 
     335               zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     336            END_3D 
    402337         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 
     338         DO_2D_00_00 
     339            zwi(ji,jj,1) = 0._wp 
     340            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
     341         END_2D 
    409342      ENDIF 
    410343      ! 
     
    416349      ! 
    417350      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 
     351         DO_2D_00_00 
     352            ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points) 
     353            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point 
     354            zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va            
     355         END_2D 
    425356         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 
     357            DO_2D_00_00 
     358               ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points) 
     359               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point 
     360               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va 
     361            END_2D 
    433362         ENDIF 
    434363      ENDIF 
     
    449378      !----------------------------------------------------------------------- 
    450379      ! 
    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 
     380      DO_3D_00_00( 2, jpkm1 ) 
     381         zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
     382      END_3D 
     383      ! 
     384      DO_2D_00_00 
     385         ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa)  
     386         pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
     387            &                                      / ( ze3va * rho0 ) * vmask(ji,jj,1)  
     388      END_2D 
     389      DO_3D_00_00( 2, jpkm1 ) 
     390         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) 
     391      END_3D 
     392      ! 
     393      DO_2D_00_00 
     394         pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1) 
     395      END_2D 
     396      DO_3DS_00_00( jpk-2, 1, -1 ) 
     397         pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - zws(ji,jj,jk) * pvv(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk) 
     398      END_3D 
    486399      ! 
    487400      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 ) 
     401         ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) / rDt - ztrdu(:,:,:) 
     402         ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) / rDt - ztrdv(:,:,:) 
     403         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt, Kmm ) 
    491404         DEALLOCATE( ztrdu, ztrdv )  
    492405      ENDIF 
    493406      !                                          ! 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' ) 
     407      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' zdf  - Ua: ', mask1=umask,               & 
     408         &                                  tab3d_2=pvv(:,:,:,Kaa), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    496409         ! 
    497410      IF( ln_timing )   CALL timing_stop('dyn_zdf') 
Note: See TracChangeset for help on using the changeset viewer.