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 13469 for NEMO/branches/2020/temporary_r4_trunk/src/OCE/DYN/dynzdf.F90 – NEMO

Ignore:
Timestamp:
2020-09-15T12:49:18+02:00 (4 years ago)
Author:
smasson
Message:

r4_trunk: first change of DO loops for routines to be merged, see #2523

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/temporary_r4_trunk/src/OCE/DYN/dynzdf.F90

    r13466 r13469  
    110110      IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity 
    111111         DO jk = 1, jpkm1 
    112             ua(:,:,jk) = ( ub(:,:,jk) + r2dt * ua(:,:,jk) ) * umask(:,:,jk) 
    113             va(:,:,jk) = ( vb(:,:,jk) + r2dt * va(:,:,jk) ) * vmask(:,:,jk) 
     112            ua(:,:,jk) = ( uu(:,:,jk,Nnn) + r2dt * ua(:,:,jk) ) * umask(:,:,jk) 
     113            va(:,:,jk) = ( vv(:,:,jk,Nnn) + r2dt * va(:,:,jk) ) * vmask(:,:,jk) 
    114114         END DO 
    115115      ELSE                                      ! applied on thickness weighted velocity 
    116116         DO jk = 1, jpkm1 
    117             ua(:,:,jk) = (         e3u_b(:,:,jk) * ub(:,:,jk)  & 
     117            ua(:,:,jk) = (         e3u_b(:,:,jk) * uu(:,:,jk,Nnn)  & 
    118118               &          + r2dt * e3u_n(:,:,jk) * ua(:,:,jk)  ) / e3u_a(:,:,jk) * umask(:,:,jk) 
    119             va(:,:,jk) = (         e3v_b(:,:,jk) * vb(:,:,jk)  & 
     119            va(:,:,jk) = (         e3v_b(:,:,jk) * vv(:,:,jk,Nnn)  & 
    120120               &          + r2dt * e3v_n(:,:,jk) * va(:,:,jk)  ) / e3v_a(:,:,jk) * vmask(:,:,jk) 
    121121         END DO 
     
    131131            va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk) 
    132132         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) 
     133         DO_2D_00_00 
     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_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) 
     137            ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) 
     138            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 
     139            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 
     140         END_2D 
     141         IF( ln_isfcav.OR.ln_drgice_imp ) THEN    ! Ocean cavities (ISF) 
     142            DO_2D_00_00 
     143               iku = miku(ji,jj)         ! top ocean level at u- and v-points  
     144               ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points) 
    137145               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) 
    138146               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 
    143          IF( ln_isfcav.OR.ln_drgice_imp ) 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,jj+1)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va 
    152                END DO 
    153             END DO 
     147               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 
     148               va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va 
     149            END_2D 
    154150         END IF 
    155151      ENDIF 
     
    162158         SELECT CASE( nldf_dyn ) 
    163159         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 
     160            DO_3D_00_00( 1, jpkm1 ) 
     161               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point 
     162               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
     163                  &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
     164               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
     165                  &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
     166               zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua 
     167               zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua 
     168               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp )  
     169               zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp ) 
     170               zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 
     171            END_3D 
    180172         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 
     173            DO_3D_00_00( 1, jpkm1 ) 
     174               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point 
     175               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
     176               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) 
     177               zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua 
     178               zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua 
     179               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 
     180               zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp ) 
     181               zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 
     182            END_3D 
    195183         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 
     184         DO_2D_00_00 
     185            zwi(ji,jj,1) = 0._wp 
     186            ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 
     187            zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) ) / ( ze3ua * e3uw_n(ji,jj,2) ) * wumask(ji,jj,2) 
     188            zWus = ( wi(ji  ,jj,2) +  wi(ji+1,jj,2) ) / ze3ua 
     189            zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp ) 
     190            zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ) ) 
     191         END_2D 
    206192      ELSE 
    207193         SELECT CASE( nldf_dyn ) 
    208194         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 
     195            DO_3D_00_00( 1, jpkm1 ) 
     196               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point 
     197               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
     198                  &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
     199               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
     200                  &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
     201               zwi(ji,jj,jk) = zzwi 
     202               zws(ji,jj,jk) = zzws 
     203               zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     204            END_3D 
    223205         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 
     206            DO_3D_00_00( 1, jpkm1 ) 
     207               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point 
     208               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
     209               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) 
     210               zwi(ji,jj,jk) = zzwi 
     211               zws(ji,jj,jk) = zzws 
     212               zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     213            END_3D 
    236214         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 
     215         DO_2D_00_00 
     216            zwi(ji,jj,1) = 0._wp 
     217            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
     218         END_2D 
    243219      ENDIF 
    244220      ! 
     
    251227      ! 
    252228      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 
     229         DO_2D_00_00 
     230            iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points 
     231            ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point 
     232            zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 
     233         END_2D 
     234         IF ( ln_isfcav.OR.ln_drgice_imp ) THEN   ! top friction (always implicit) 
     235            DO_2D_00_00 
     236               !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed 
     237               iku = miku(ji,jj)       ! ocean top level at u- and v-points  
    256238               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 
    260          IF ( ln_isfcav.OR.ln_drgice_imp ) 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 
     239               zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 
     240            END_2D 
    269241         END IF 
    270242      ENDIF 
     
    285257      !----------------------------------------------------------------------- 
    286258      ! 
    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 
     259      DO_3D_00_00( 2, jpkm1 ) 
     260         zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
     261      END_3D 
     262      ! 
     263      DO_2D_00_00 
     264         ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1)  
     265         ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
     266            &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1)  
     267      END_2D 
     268      DO_3D_00_00( 2, jpkm1 ) 
     269         ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 
     270      END_3D 
     271      ! 
     272      DO_2D_00_00 
     273         ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
     274      END_2D 
     275      DO_3D_00_00( jpk-2, 1, -1 ) 
     276         ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
     277      END_3D 
    322278      ! 
    323279      !              !==  Vertical diffusion on v  ==! 
     
    328284         SELECT CASE( nldf_dyn ) 
    329285         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 
     286            DO_3D_00_00( 1, jpkm1 ) 
     287               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point 
     288               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
     289                  &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
     290               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
     291                  &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
     292               zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va 
     293               zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va 
     294               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp ) 
     295               zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp ) 
     296               zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 
     297            END_3D 
    346298         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 
     299            DO_3D_00_00( 1, jpkm1 ) 
     300               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point 
     301               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
     302               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) 
     303               zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va 
     304               zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va 
     305               zwi(ji,jj,jk) = zzwi  + zdt * MIN( zWvi, 0._wp ) 
     306               zws(ji,jj,jk) = zzws  - zdt * MAX( zWvs, 0._wp ) 
     307               zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 
     308            END_3D 
    361309         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 
     310         DO_2D_00_00 
     311            zwi(ji,jj,1) = 0._wp 
     312            ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 
     313            zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw_n(ji,jj,2) ) * wvmask(ji,jj,2) 
     314            zWvs = ( wi(ji,jj  ,2) +  wi(ji,jj+1,2) ) / ze3va 
     315            zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp ) 
     316            zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) ) 
     317         END_2D 
    372318      ELSE 
    373319         SELECT CASE( nldf_dyn ) 
    374320         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 
     321            DO_3D_00_00( 1, jpkm1 ) 
     322               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point 
     323               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
     324                  &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
     325               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
     326                  &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
     327               zwi(ji,jj,jk) = zzwi 
     328               zws(ji,jj,jk) = zzws 
     329               zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     330            END_3D 
    389331         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 
     332            DO_3D_00_00( 1, jpkm1 ) 
     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  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
     335               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) 
     336               zwi(ji,jj,jk) = zzwi 
     337               zws(ji,jj,jk) = zzws 
     338               zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     339            END_3D 
    402340         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 
     341         DO_2D_00_00 
     342            zwi(ji,jj,1) = 0._wp 
     343            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
     344         END_2D 
    409345      ENDIF 
    410346      ! 
     
    416352      ! 
    417353      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) 
     354         DO_2D_00_00 
     355            ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points) 
     356            ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point 
     357            zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va            
     358         END_2D 
     359         IF ( ln_isfcav.OR.ln_drgice_imp ) THEN 
     360            DO_2D_00_00 
     361               ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points) 
    421362               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 
    425          IF ( ln_isfcav.OR.ln_drgice_imp ) 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,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va 
    431                END DO 
    432             END DO 
     363               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va 
     364            END_2D 
    433365         ENDIF 
    434366      ENDIF 
     
    449381      !----------------------------------------------------------------------- 
    450382      ! 
    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 
     383      DO_3D_00_00( 2, jpkm1 ) 
     384         zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
     385      END_3D 
     386      ! 
     387      DO_2D_00_00 
     388         ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1)  
     389         va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
     390            &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1)  
     391      END_2D 
     392      DO_3D_00_00( 2, jpkm1 ) 
     393         va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 
     394      END_3D 
     395      ! 
     396      DO_2D_00_00 
     397         va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
     398      END_2D 
     399      DO_3D_00_00( jpk-2, 1, -1 ) 
     400         va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
     401      END_3D 
    486402      ! 
    487403      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics 
    488          ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:) 
    489          ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:) 
     404         ztrdu(:,:,:) = ( ua(:,:,:) - uu(:,:,:,Nnn) ) / r2dt - ztrdu(:,:,:) 
     405         ztrdv(:,:,:) = ( va(:,:,:) - vv(:,:,:,Nnn) ) / r2dt - ztrdv(:,:,:) 
    490406         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt ) 
    491407         DEALLOCATE( ztrdu, ztrdv )  
Note: See TracChangeset for help on using the changeset viewer.