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 10180 for NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_adv_umx.F90 – NEMO

Ignore:
Timestamp:
2018-10-08T17:06:04+02:00 (5 years ago)
Author:
smasson
Message:

dev_r10164_HPC09_ESIWACE_PREP_MERGE: action 4b: reduce communications in si3, see #2133

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_adv_umx.F90

    r10170 r10180  
    7171      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices 
    7272      INTEGER  ::   initad                  ! number of sub-timestep for the advection 
     73      INTEGER  ::   ipl                     ! third dimention of tracer array 
     74 
    7375      REAL(wp) ::   zcfl , zusnit, zdt      !   -      - 
    74       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zudy, zvdx, zcu_box, zcv_box 
     76      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zudy, zvdx, zcu_box, zcv_box 
     77      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zpato 
    7578      !!---------------------------------------------------------------------- 
    7679      ! 
     
    7881      ! 
    7982      ALLOCATE( zudy(jpi,jpj) , zvdx(jpi,jpj) , zcu_box(jpi,jpj) , zcv_box(jpi,jpj) ) 
     83      ALLOCATE( zpato(jpi,jpj,1) ) 
    8084      ! 
    8185      ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !         
     
    109113      END DO 
    110114 
     115      zpato (:,:,1) = pato_i(:,:) 
     116 
    111117      !---------------! 
    112118      !== advection ==! 
    113119      !---------------! 
    114120      DO jt = 1, initad 
    115          CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:) )             ! Open water area  
    116          DO jl = 1, jpl 
    117             CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl) )         ! Ice area 
    118             CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_i(:,:,jl) )         ! Ice  volume 
    119             CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, psv_i(:,:,jl) )        ! Salt content 
    120             CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, poa_i(:,:,jl) )        ! Age content 
    121             DO jk = 1, nlay_i 
    122                CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_i(:,:,jk,jl) )   ! Ice  heat content 
    123             END DO 
    124             CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,jl) )         ! Snow volume 
    125             DO jk = 1, nlay_s 
    126                CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,jk,jl) )   ! Snow heat content 
    127             END DO 
    128             IF ( ln_pnd_H12 ) THEN 
    129                CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl) )     ! Melt pond fraction 
    130                CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,jl) )     ! Melt pond volume 
    131             ENDIF 
    132          END DO 
    133       END DO 
    134       ! 
    135       DEALLOCATE( zudy, zvdx, zcu_box, zcv_box ) 
     121         CALL adv_umx( k_order, kt,   1, zdt, zudy, zvdx, zcu_box, zcv_box, zpato(:,:,1) )        ! Open water area  
     122         CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,:) )         ! Ice area 
     123         CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pv_i(:,:,:) )         ! Ice  volume 
     124         CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, psv_i(:,:,:) )        ! Salt content 
     125         CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, poa_i(:,:,:) )        ! Age content 
     126         DO jk = 1, nlay_i 
     127            CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pe_i(:,:,jk,:) )   ! Ice  heat content 
     128         END DO 
     129         CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,:) )         ! Snow volume 
     130         DO jk = 1, nlay_s 
     131            CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,jk,:) )   ! Snow heat content 
     132         END DO 
     133         IF ( ln_pnd_H12 ) THEN 
     134            CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,:) )     ! Melt pond fraction 
     135            CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,:) )     ! Melt pond volume 
     136         ENDIF 
     137      END DO 
     138      ! 
     139      pato_i(:,:) = zpato (:,:,1) 
     140      ! 
     141      DEALLOCATE( zudy, zvdx, zcu_box, zcv_box, zpato ) 
    136142      ! 
    137143   END SUBROUTINE ice_dyn_adv_umx 
    138144 
    139145    
    140    SUBROUTINE adv_umx( k_order, kt, pdt, puc, pvc, pubox, pvbox, ptc ) 
     146   SUBROUTINE adv_umx( k_order, kt, ipl, pdt, puc, pvc, pubox, pvbox, ptc ) 
    141147      !!---------------------------------------------------------------------- 
    142148      !!                  ***  ROUTINE adv_umx  *** 
     
    151157      !! ** Action : - pt  the after advective tracer 
    152158      !!---------------------------------------------------------------------- 
    153       INTEGER                     , INTENT(in   ) ::   k_order        ! order of the ULTIMATE scheme 
    154       INTEGER                     , INTENT(in   ) ::   kt             ! number of iteration 
    155       REAL(wp)                    , INTENT(in   ) ::   pdt            ! tracer time-step 
    156       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc  , pvc     ! 2 ice velocity components => u*e2 
    157       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pubox, pvbox   ! upstream velocity 
    158       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ptc            ! tracer content field 
    159       ! 
    160       INTEGER  ::   ji, jj           ! dummy loop indices   
     159      INTEGER                         , INTENT(in   ) ::   k_order        ! order of the ULTIMATE scheme 
     160      INTEGER                         , INTENT(in   ) ::   kt             ! number of iteration 
     161      INTEGER                         , INTENT(in   ) ::   ipl            ! third dimension of tracer array 
     162      REAL(wp)                        , INTENT(in   ) ::   pdt            ! tracer time-step 
     163      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   puc  , pvc     ! 2 ice velocity components => u*e2 
     164      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   pubox, pvbox   ! upstream velocity 
     165      REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(inout) ::   ptc            ! tracer content field 
     166      ! 
     167      INTEGER  ::   ji, jj, jl       ! dummy loop indices   
     168 
    161169      REAL(wp) ::   ztra             ! local scalar 
    162       REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ups, zfu_ho, zt_u, zt_ups 
    163       REAL(wp), DIMENSION(jpi,jpj) ::   zfv_ups, zfv_ho, zt_v, ztrd 
     170      REAL(wp), DIMENSION(jpi,jpj,ipl) ::   zfu_ups, zfu_ho, zt_u, zt_ups 
     171      REAL(wp), DIMENSION(jpi,jpj,ipl) ::   zfv_ups, zfv_ho, zt_v, ztrd 
     172 
     173      DO jl = 1, ipl 
    164174      !!---------------------------------------------------------------------- 
    165175      ! 
     
    168178      DO jj = 1, jpjm1         ! upstream tracer flux in the i and j direction 
    169179         DO ji = 1, fs_jpim1   ! vector opt. 
    170             zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * ptc(ji,jj) + MIN( puc(ji,jj), 0._wp ) * ptc(ji+1,jj) 
    171             zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * ptc(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * ptc(ji,jj+1) 
     180            zfu_ups(ji,jj,jl) = MAX( puc(ji,jj), 0._wp ) * ptc(ji,jj,jl) + MIN( puc(ji,jj), 0._wp ) * ptc(ji+1,jj,jl) 
     181            zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj), 0._wp ) * ptc(ji,jj,jl) + MIN( pvc(ji,jj), 0._wp ) * ptc(ji,jj+1,jl) 
    172182         END DO 
    173183      END DO 
    174184       
    175       DO jj = 2, jpjm1            ! total intermediate advective trends 
    176          DO ji = fs_2, fs_jpim1   ! vector opt. 
    177             ztra = - (   zfu_ups(ji,jj) - zfu_ups(ji-1,jj  )   & 
    178                &       + zfv_ups(ji,jj) - zfv_ups(ji  ,jj-1)   ) * r1_e1e2t(ji,jj) 
     185         DO jj = 2, jpjm1            ! total intermediate advective trends 
     186            DO ji = fs_2, fs_jpim1   ! vector opt. 
     187               ztra = - (   zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj  ,jl)   & 
     188                  &       + zfv_ups(ji,jj,jl) - zfv_ups(ji  ,jj-1,jl)   ) * r1_e1e2t(ji,jj) 
    179189            ! 
    180             ztrd(ji,jj) =                         ztra                         ! upstream trend [ -div(uh) or -div(uhT) ]   
    181             zt_ups (ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)      ! guess after content field with monotonic scheme 
    182          END DO 
    183       END DO 
    184       CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign) 
     190               ztrd(ji,jj,jl) =                         ztra                         ! upstream trend [ -div(uh) or -div(uhT) ]   
     191               zt_ups (ji,jj,jl) = ( ptc(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1)   ! guess after content field with monotonic scheme 
     192            END DO 
     193         END DO 
     194      END DO 
    185195       
    186196      ! High order (_ho) fluxes  
     
    188198      SELECT CASE( k_order ) 
    189199      CASE ( 20 )                          ! centered second order 
    190          DO jj = 1, jpjm1 
    191             DO ji = 1, fs_jpim1   ! vector opt. 
    192                zfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( ptc(ji,jj) + ptc(ji+1,jj) ) 
    193                zfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( ptc(ji,jj) + ptc(ji,jj+1) ) 
     200         DO jl = 1, ipl 
     201            DO jj = 1, jpjm1 
     202               DO ji = 1, fs_jpim1   ! vector opt. 
     203                  zfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj) * ( ptc(ji,jj,jl) + ptc(ji+1,jj,jl) ) 
     204                  zfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj) * ( ptc(ji,jj,jl) + ptc(ji,jj+1,jl) ) 
     205               END DO 
    194206            END DO 
    195207         END DO 
    196208         ! 
    197209      CASE ( 1:5 )                      ! 1st to 5th order ULTIMATE-MACHO scheme 
    198          CALL macho( k_order, kt, pdt, ptc, puc, pvc, pubox, pvbox, zt_u, zt_v ) 
    199          ! 
    200          DO jj = 1, jpjm1 
    201             DO ji = 1, fs_jpim1   ! vector opt. 
    202                zfu_ho(ji,jj) = puc(ji,jj) * zt_u(ji,jj) 
    203                zfv_ho(ji,jj) = pvc(ji,jj) * zt_v(ji,jj) 
     210         CALL macho( k_order, kt, ipl, pdt, ptc, puc, pvc, pubox, pvbox, zt_u, zt_v ) 
     211         ! 
     212         DO jl = 1, ipl 
     213            DO jj = 1, jpjm1 
     214               DO ji = 1, fs_jpim1   ! vector opt. 
     215                  zfu_ho(ji,jj,jl) = puc(ji,jj) * zt_u(ji,jj,jl) 
     216                  zfv_ho(ji,jj,jl) = pvc(ji,jj) * zt_v(ji,jj,jl) 
     217               END DO 
    204218            END DO 
    205219         END DO 
     
    209223      ! antidiffusive flux : high order minus low order 
    210224      ! -------------------------------------------------- 
    211       DO jj = 1, jpjm1 
    212          DO ji = 1, fs_jpim1   ! vector opt. 
    213             zfu_ho(ji,jj) = zfu_ho(ji,jj) - zfu_ups(ji,jj) 
    214             zfv_ho(ji,jj) = zfv_ho(ji,jj) - zfv_ups(ji,jj) 
    215          END DO 
    216       END DO 
     225      DO jl = 1, ipl 
     226         DO jj = 1, jpjm1 
     227            DO ji = 1, fs_jpim1   ! vector opt. 
     228               zfu_ho(ji,jj,jl) = zfu_ho(ji,jj,jl) - zfu_ups(ji,jj,jl) 
     229               zfv_ho(ji,jj,jl) = zfv_ho(ji,jj,jl) - zfv_ups(ji,jj,jl) 
     230            END DO 
     231         END DO 
     232      END DO 
     233 
     234      CALL lbc_lnk("icedyn_adv_umx",zt_ups, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign) 
    217235       
    218236      ! monotonicity algorithm 
    219237      ! ------------------------- 
    220       CALL nonosc_2d( ptc, zfu_ho, zfv_ho, zt_ups, pdt ) 
     238      CALL nonosc_2d( ipl, ptc, zfu_ho, zfv_ho, zt_ups, pdt ) 
    221239       
    222240      ! final trend with corrected fluxes 
    223241      ! ------------------------------------ 
    224       DO jj = 2, jpjm1 
    225          DO ji = fs_2, fs_jpim1   ! vector opt.   
    226             ztra       = ztrd(ji,jj)  - (  zfu_ho(ji,jj) - zfu_ho(ji-1,jj  )   & 
    227                &                         + zfv_ho(ji,jj) - zfv_ho(ji  ,jj-1) ) * r1_e1e2t(ji,jj)   
    228             ptc(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 
     242      DO jl = 1, ipl 
     243         DO jj = 2, jpjm1 
     244            DO ji = fs_2, fs_jpim1   ! vector opt.   
     245               ztra       = ztrd(ji,jj,jl)  - (  zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj  ,jl)   & 
     246                  &                         +    zfv_ho(ji,jj,jl) - zfv_ho(ji  ,jj-1,jl) ) * r1_e1e2t(ji,jj)   
     247               ptc(ji,jj,jl) = ptc(ji,jj,jl) + pdt * ztra * tmask(ji,jj,1) 
     248            END DO 
    229249         END DO 
    230250      END DO 
     
    233253   END SUBROUTINE adv_umx 
    234254 
    235  
    236    SUBROUTINE macho( k_order, kt, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v ) 
     255   SUBROUTINE macho( k_order, kt, ipl, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v ) 
    237256      !!--------------------------------------------------------------------- 
    238257      !!                    ***  ROUTINE ultimate_x  *** 
     
    245264      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    246265      !!---------------------------------------------------------------------- 
    247       INTEGER                     , INTENT(in   ) ::   k_order    ! order of the ULTIMATE scheme 
    248       INTEGER                     , INTENT(in   ) ::   kt         ! number of iteration 
    249       REAL(wp)                    , INTENT(in   ) ::   pdt        ! tracer time-step 
    250       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptc        ! tracer fields 
    251       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc, pvc   ! 2 ice velocity components 
    252       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pubox, pvbox   ! upstream velocity 
    253       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u, pt_v ! tracer at u- and v-points  
    254       ! 
    255       INTEGER  ::   ji, jj    ! dummy loop indices 
     266      INTEGER                         , INTENT(in   ) ::   k_order    ! order of the ULTIMATE scheme 
     267      INTEGER                         , INTENT(in   ) ::   kt         ! number of iteration 
     268      INTEGER                         , INTENT(in   ) ::   ipl        ! third dimension of tracer array 
     269      REAL(wp)                        , INTENT(in   ) ::   pdt        ! tracer time-step 
     270      REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(in   ) ::   ptc        ! tracer fields 
     271      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   puc, pvc   ! 2 ice velocity components 
     272      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   pubox, pvbox   ! upstream velocity 
     273      REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(  out) ::   pt_u, pt_v ! tracer at u- and v-points  
     274      ! 
     275      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
    256276      REAL(wp) ::   zc_box    !   -      - 
    257       REAL(wp), DIMENSION(jpi,jpj) :: zzt 
     277      REAL(wp), DIMENSION(jpi,jpj,ipl) :: zzt 
    258278      !!---------------------------------------------------------------------- 
    259279      ! 
     
    261281         ! 
    262282         !                                                           !--  ultimate interpolation of pt at u-point  --! 
    263          CALL ultimate_x( k_order, pdt, ptc, puc, pt_u ) 
     283         CALL ultimate_x( k_order, ipl, pdt, ptc, puc, pt_u ) 
    264284         ! 
    265285         !                                                           !--  advective form update in zzt  --! 
    266          DO jj = 2, jpjm1 
    267             DO ji = fs_2, fs_jpim1   ! vector opt. 
    268                zzt(ji,jj) = ptc(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj)  & 
    269                   &                    - ptc  (ji,jj) * pdt * ( puc (ji,jj) - puc (ji-1,jj) ) * r1_e1e2t(ji,jj) 
    270                zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 
     286         DO jl = 1, ipl 
     287            DO jj = 2, jpjm1 
     288               DO ji = fs_2, fs_jpim1   ! vector opt. 
     289                  zzt(ji,jj,jl) = ptc(ji,jj,jl) - pubox(ji,jj) * pdt * ( pt_u(ji,jj,jl) - pt_u(ji-1,jj,jl) ) * r1_e1t(ji,jj)  & 
     290                   &            - ptc(ji,jj,jl) * pdt * ( puc (ji,jj) - puc (ji-1,jj) ) * r1_e1e2t(ji,jj) 
     291                  zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1) 
     292               END DO 
    271293            END DO 
    272294         END DO 
     
    274296         ! 
    275297         !                                                           !--  ultimate interpolation of pt at v-point  --! 
    276          CALL ultimate_y( k_order, pdt, zzt, pvc, pt_v ) 
     298         CALL ultimate_y( k_order, ipl, pdt, zzt, pvc, pt_v ) 
    277299         ! 
    278300      ELSE                                                  !==  even ice time step:  adv_y then adv_x  ==! 
    279301         ! 
    280302         !                                                           !--  ultimate interpolation of pt at v-point  --! 
    281          CALL ultimate_y( k_order, pdt, ptc, pvc, pt_v ) 
     303         CALL ultimate_y( k_order, ipl, pdt, ptc, pvc, pt_v ) 
    282304         ! 
    283305         !                                                           !--  advective form update in zzt  --! 
    284          DO jj = 2, jpjm1 
    285             DO ji = fs_2, fs_jpim1 
    286                zzt(ji,jj) = ptc(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj)  & 
    287                   &                    - ptc  (ji,jj) * pdt * ( pvc (ji,jj) - pvc (ji,jj-1) ) * r1_e1e2t(ji,jj) 
    288                zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 
     306         DO jl = 1, ipl 
     307            DO jj = 2, jpjm1 
     308               DO ji = fs_2, fs_jpim1 
     309                  zzt(ji,jj,jl) = ptc(ji,jj,jl) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj,jl) - pt_v(ji,jj-1,jl) ) * r1_e2t(ji,jj)  & 
     310                     &                    - ptc  (ji,jj,jl) * pdt * ( pvc (ji,jj) - pvc (ji,jj-1) ) * r1_e1e2t(ji,jj) 
     311                  zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1) 
     312               END DO 
    289313            END DO 
    290314         END DO 
     
    292316         ! 
    293317         !                                                           !--  ultimate interpolation of pt at u-point  --! 
    294          CALL ultimate_x( k_order, pdt, zzt, puc, pt_u ) 
     318         CALL ultimate_x( k_order, ipl, pdt, zzt, puc, pt_u ) 
    295319         !       
    296320      ENDIF       
     
    299323 
    300324 
    301    SUBROUTINE ultimate_x( k_order, pdt, pt, puc, pt_u ) 
     325   SUBROUTINE ultimate_x( k_order, ipl, pdt, pt, puc, pt_u ) 
    302326      !!--------------------------------------------------------------------- 
    303327      !!                    ***  ROUTINE ultimate_x  *** 
     
    310334      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    311335      !!---------------------------------------------------------------------- 
    312       INTEGER                     , INTENT(in   ) ::   k_order   ! ocean time-step index 
    313       REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step 
    314       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc       ! ice i-velocity component 
    315       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt        ! tracer fields 
    316       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u      ! tracer at u-point  
    317       ! 
    318       INTEGER  ::   ji, jj       ! dummy loop indices 
     336      INTEGER                         , INTENT(in   ) ::   k_order   ! ocean time-step index 
     337      INTEGER                         , INTENT(in   ) ::   ipl        ! third dimension of tracer array 
     338      REAL(wp)                        , INTENT(in   ) ::   pdt       ! tracer time-step 
     339      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   puc       ! ice i-velocity component 
     340      REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(in   ) ::   pt        ! tracer fields 
     341      REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(  out) ::   pt_u      ! tracer at u-point  
     342      ! 
     343      INTEGER  ::   ji, jj, jl       ! dummy loop indices 
    319344      REAL(wp) ::   zcu, zdx2, zdx4    !   -      - 
    320       REAL(wp), DIMENSION(jpi,jpj) :: ztu1, ztu2, ztu3, ztu4 
     345      REAL(wp), DIMENSION(jpi,jpj) :: ztu1, ztu3 
     346      REAL(wp), DIMENSION(jpi,jpj,ipl) :: ztu2, ztu4 
    321347      !!---------------------------------------------------------------------- 
    322348      ! 
    323349      !                                                     !--  Laplacian in i-direction  --! 
    324       DO jj = 2, jpjm1         ! First derivative (gradient) 
    325          DO ji = 1, fs_jpim1 
    326             ztu1(ji,jj) = ( pt(ji+1,jj) - pt(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 
    327          END DO 
    328          !                     ! Second derivative (Laplacian) 
    329          DO ji = fs_2, fs_jpim1 
    330             ztu2(ji,jj) = ( ztu1(ji,jj) - ztu1(ji-1,jj) ) * r1_e1t(ji,jj) 
     350      DO jl = 1, ipl 
     351         DO jj = 2, jpjm1         ! First derivative (gradient) 
     352            DO ji = 1, fs_jpim1 
     353               ztu1(ji,jj) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 
     354            END DO 
     355            !                     ! Second derivative (Laplacian) 
     356            DO ji = fs_2, fs_jpim1 
     357               ztu2(ji,jj,jl) = ( ztu1(ji,jj) - ztu1(ji-1,jj) ) * r1_e1t(ji,jj) 
     358            END DO 
    331359         END DO 
    332360      END DO 
     
    334362      ! 
    335363      !                                                     !--  BiLaplacian in i-direction  --! 
    336       DO jj = 2, jpjm1         ! Third derivative 
    337          DO ji = 1, fs_jpim1 
    338             ztu3(ji,jj) = ( ztu2(ji+1,jj) - ztu2(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 
    339          END DO 
    340          !                     ! Fourth derivative 
    341          DO ji = fs_2, fs_jpim1 
    342             ztu4(ji,jj) = ( ztu3(ji,jj) - ztu3(ji-1,jj) ) * r1_e1t(ji,jj) 
     364      DO jl = 1, ipl 
     365         DO jj = 2, jpjm1         ! Third derivative 
     366            DO ji = 1, fs_jpim1 
     367               ztu3(ji,jj) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 
     368            END DO 
     369            !                     ! Fourth derivative 
     370            DO ji = fs_2, fs_jpim1 
     371               ztu4(ji,jj,jl) = ( ztu3(ji,jj) - ztu3(ji-1,jj) ) * r1_e1t(ji,jj) 
     372            END DO 
    343373         END DO 
    344374      END DO 
     
    350380      CASE( 1 )                                                   !==  1st order central TIM  ==! (Eq. 21) 
    351381         !         
    352          DO jj = 2, jpjm1 
    353             DO ji = 1, fs_jpim1   ! vector opt. 
    354                pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                               pt(ji+1,jj) + pt(ji,jj)   & 
    355                   &                                    - SIGN( 1._wp, puc(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) ) 
     382         DO jl = 1, ipl 
     383            DO jj = 2, jpjm1 
     384               DO ji = 1, fs_jpim1   ! vector opt. 
     385                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                               pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
     386                     &                                       - SIGN( 1._wp, puc(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
     387               END DO 
    356388            END DO 
    357389         END DO 
     
    359391      CASE( 2 )                                                   !==  2nd order central TIM  ==! (Eq. 23) 
    360392         ! 
    361          DO jj = 2, jpjm1 
    362             DO ji = 1, fs_jpim1   ! vector opt. 
    363                zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    364                pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                                   pt(ji+1,jj) + pt(ji,jj)   & 
    365                   &                                               -              zcu   * ( pt(ji+1,jj) - pt(ji,jj) ) )  
     393         DO jl = 1, ipl 
     394            DO jj = 2, jpjm1 
     395               DO ji = 1, fs_jpim1   ! vector opt. 
     396                  zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     397                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                               pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
     398                     &                                              -              zcu   * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )  
     399               END DO 
    366400            END DO 
    367401         END DO 
     
    369403      CASE( 3 )                                                   !==  3rd order central TIM  ==! (Eq. 24) 
    370404         ! 
    371          DO jj = 2, jpjm1 
    372             DO ji = 1, fs_jpim1   ! vector opt. 
    373                zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    374                zdx2 = e1u(ji,jj) * e1u(ji,jj) 
     405         DO jl = 1, ipl 
     406            DO jj = 2, jpjm1 
     407               DO ji = 1, fs_jpim1   ! vector opt. 
     408                  zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     409                  zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    375410!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj) 
    376                pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (         (                         pt  (ji+1,jj) + pt  (ji,jj)        & 
    377                   &                                               -              zcu   * ( pt  (ji+1,jj) - pt  (ji,jj) )  )   & 
    378                   &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                         ztu2(ji+1,jj) + ztu2(ji,jj)        & 
    379                   &                                               - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj) - ztu2(ji,jj) )  )   ) 
     411                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                     pt  (ji+1,jj,jl) + pt  (ji,jj,jl)        & 
     412                     &                                              -              zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) )  )   & 
     413                     &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                        ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)        & 
     414                     &                                              - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) )  )   ) 
     415               END DO 
    380416            END DO 
    381417         END DO 
     
    383419      CASE( 4 )                                                   !==  4th order central TIM  ==! (Eq. 27) 
    384420         ! 
    385          DO jj = 2, jpjm1 
    386             DO ji = 1, fs_jpim1   ! vector opt. 
    387                zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    388                zdx2 = e1u(ji,jj) * e1u(ji,jj) 
     421         DO jl = 1, ipl 
     422            DO jj = 2, jpjm1 
     423               DO ji = 1, fs_jpim1   ! vector opt. 
     424                  zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     425                  zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    389426!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj) 
    390                pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (         (                   pt  (ji+1,jj) + pt  (ji,jj)        & 
    391                   &                                               -          zcu * ( pt  (ji+1,jj) - pt  (ji,jj) )  )   & 
    392                   &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                   ztu2(ji+1,jj) + ztu2(ji,jj)        & 
    393                   &                                               - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) )  )   ) 
     427                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                     pt  (ji+1,jj,jl) + pt  (ji,jj,jl)        & 
     428                     &                                                    -          zcu * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) )  )   & 
     429                     &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                        ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)        & 
     430                     &                                                    - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) )  )   ) 
     431               END DO 
    394432            END DO 
    395433         END DO 
     
    397435      CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29) 
    398436         ! 
    399          DO jj = 2, jpjm1 
    400             DO ji = 1, fs_jpim1   ! vector opt. 
    401                zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    402                zdx2 = e1u(ji,jj) * e1u(ji,jj) 
     437         DO jl = 1, ipl 
     438            DO jj = 2, jpjm1 
     439               DO ji = 1, fs_jpim1   ! vector opt. 
     440                  zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     441                  zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    403442!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj) 
    404                zdx4 = zdx2 * zdx2 
    405                pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (               (                   pt  (ji+1,jj) + pt  (ji,jj)       & 
    406                   &                                                     -          zcu * ( pt  (ji+1,jj) - pt  (ji,jj) ) )   & 
    407                   &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) *     (                   ztu2(ji+1,jj) + ztu2(ji,jj)       & 
    408                   &                                                     - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) ) )   & 
    409                   &        + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj) + ztu4(ji,jj)       & 
    410                   &                                               - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj) - ztu4(ji,jj) ) ) ) 
     443                  zdx4 = zdx2 * zdx2 
     444                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (           (                   pt  (ji+1,jj,jl) + pt  (ji,jj,jl)       & 
     445                     &                                                    -          zcu * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) )   & 
     446                     &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) *    (                   ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)       & 
     447                     &                                                    - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) )   & 
     448                     &       + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl)       & 
     449                     &                                              - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) ) 
     450               END DO 
    411451            END DO 
    412452         END DO 
     
    417457    
    418458  
    419    SUBROUTINE ultimate_y( k_order, pdt, pt, pvc, pt_v ) 
     459   SUBROUTINE ultimate_y( k_order, ipl, pdt, pt, pvc, pt_v ) 
    420460      !!--------------------------------------------------------------------- 
    421461      !!                    ***  ROUTINE ultimate_y  *** 
     
    428468      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    429469      !!---------------------------------------------------------------------- 
    430       INTEGER                     , INTENT(in   ) ::   k_order   ! ocean time-step index 
    431       REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step 
    432       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvc       ! ice j-velocity component 
    433       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt        ! tracer fields 
    434       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_v      ! tracer at v-point  
    435       ! 
    436       INTEGER  ::   ji, jj       ! dummy loop indices 
     470      INTEGER                         , INTENT(in   ) ::   k_order   ! ocean time-step index 
     471      INTEGER                         , INTENT(in   ) ::   ipl       ! third dimension of tracer array 
     472      REAL(wp)                        , INTENT(in   ) ::   pdt       ! tracer time-step 
     473      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   pvc       ! ice j-velocity component 
     474      REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(in   ) ::   pt        ! tracer fields 
     475      REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(  out) ::   pt_v      ! tracer at v-point  
     476      ! 
     477      INTEGER  ::   ji, jj, jl       ! dummy loop indices 
    437478      REAL(wp) ::   zcv, zdy2, zdy4    !   -      - 
    438       REAL(wp), DIMENSION(jpi,jpj) :: ztv1, ztv2, ztv3, ztv4 
     479      REAL(wp), DIMENSION(jpi,jpj) :: ztv1, ztv3 
     480      REAL(wp), DIMENSION(jpi,jpj,ipl) :: ztv2, ztv4 
    439481      !!---------------------------------------------------------------------- 
    440482      ! 
    441483      !                                                     !--  Laplacian in j-direction  --! 
    442       DO jj = 1, jpjm1         ! First derivative (gradient) 
    443          DO ji = fs_2, fs_jpim1 
    444             ztv1(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 
    445          END DO 
    446       END DO 
    447       DO jj = 2, jpjm1         ! Second derivative (Laplacian) 
    448          DO ji = fs_2, fs_jpim1 
    449             ztv2(ji,jj) = ( ztv1(ji,jj) - ztv1(ji,jj-1) ) * r1_e2t(ji,jj) 
     484      DO jl = 1, ipl 
     485         DO jj = 1, jpjm1         ! First derivative (gradient) 
     486            DO ji = fs_2, fs_jpim1 
     487               ztv1(ji,jj) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 
     488            END DO 
     489         END DO 
     490         DO jj = 2, jpjm1         ! Second derivative (Laplacian) 
     491            DO ji = fs_2, fs_jpim1 
     492               ztv2(ji,jj,jl) = ( ztv1(ji,jj) - ztv1(ji,jj-1) ) * r1_e2t(ji,jj) 
     493            END DO 
    450494         END DO 
    451495      END DO 
     
    453497      ! 
    454498      !                                                     !--  BiLaplacian in j-direction  --! 
    455       DO jj = 1, jpjm1         ! First derivative 
    456          DO ji = fs_2, fs_jpim1 
    457             ztv3(ji,jj) = ( ztv2(ji,jj+1) - ztv2(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 
    458          END DO 
    459       END DO 
    460       DO jj = 2, jpjm1         ! Second derivative 
    461          DO ji = fs_2, fs_jpim1 
    462             ztv4(ji,jj) = ( ztv3(ji,jj) - ztv3(ji,jj-1) ) * r1_e2t(ji,jj) 
     499      DO jl = 1, ipl 
     500         DO jj = 1, jpjm1         ! First derivative 
     501            DO ji = fs_2, fs_jpim1 
     502            ztv3(ji,jj) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 
     503            END DO 
     504         END DO 
     505         DO jj = 2, jpjm1         ! Second derivative 
     506            DO ji = fs_2, fs_jpim1 
     507               ztv4(ji,jj,jl) = ( ztv3(ji,jj) - ztv3(ji,jj-1) ) * r1_e2t(ji,jj) 
     508            END DO 
    463509         END DO 
    464510      END DO 
     
    469515      ! 
    470516      CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21) 
    471          DO jj = 1, jpjm1 
    472             DO ji = fs_2, fs_jpim1 
    473                pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1) + pt(ji,jj) )  & 
    474                   &                                     - SIGN( 1._wp, pvc(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) ) 
     517         DO jl = 1, ipl 
     518            DO jj = 1, jpjm1 
     519               DO ji = fs_2, fs_jpim1 
     520                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                           ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  & 
     521                     &                                     - SIGN( 1._wp, pvc(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
     522               END DO 
    475523            END DO 
    476524         END DO 
    477525         ! 
    478526      CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23) 
    479          DO jj = 1, jpjm1 
    480             DO ji = fs_2, fs_jpim1 
    481                zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    482                pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (        ( pt(ji,jj+1) + pt(ji,jj) )  & 
    483                   &                                     - zcv * ( pt(ji,jj+1) - pt(ji,jj) ) ) 
     527         DO jl = 1, ipl 
     528            DO jj = 1, jpjm1 
     529               DO ji = fs_2, fs_jpim1 
     530                  zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     531                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (     ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  & 
     532                     &                                     - zcv * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
     533               END DO 
    484534            END DO 
    485535         END DO 
     
    487537         ! 
    488538      CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24) 
    489          DO jj = 1, jpjm1 
    490             DO ji = fs_2, fs_jpim1 
    491                zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    492                zdy2 = e2v(ji,jj) * e2v(ji,jj) 
     539         DO jl = 1, ipl 
     540            DO jj = 1, jpjm1 
     541               DO ji = fs_2, fs_jpim1 
     542                  zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     543                  zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    493544!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj) 
    494                pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                                 ( pt  (ji,jj+1) + pt  (ji,jj)       & 
    495                   &                                     -                        zcv   * ( pt  (ji,jj+1) - pt  (ji,jj) ) )   & 
    496                   &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1) + ztv2(ji,jj)       & 
    497                   &                                               - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) ) 
     545                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)       & 
     546                     &                                     -                        zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) )   & 
     547                     &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)       & 
     548                     &                                               - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) ) 
     549               END DO 
    498550            END DO 
    499551         END DO 
    500552         ! 
    501553      CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27) 
    502          DO jj = 1, jpjm1 
    503             DO ji = fs_2, fs_jpim1 
    504                zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    505                zdy2 = e2v(ji,jj) * e2v(ji,jj) 
     554         DO jl = 1, ipl 
     555            DO jj = 1, jpjm1 
     556               DO ji = fs_2, fs_jpim1 
     557                  zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     558                  zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    506559!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj) 
    507                pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                           ( pt  (ji,jj+1) + pt  (ji,jj)       & 
    508                   &                                               -          zcv * ( pt  (ji,jj+1) - pt  (ji,jj) ) )   & 
    509                   &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                   ztv2(ji,jj+1) + ztv2(ji,jj)       & 
    510                   &                                               - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) ) 
     560                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                        ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)       & 
     561                     &                                               -          zcv * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) )   & 
     562                     &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                   ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)       & 
     563                     &                                               - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) ) 
     564               END DO 
    511565            END DO 
    512566         END DO 
    513567         ! 
    514568      CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29) 
    515          DO jj = 1, jpjm1 
    516             DO ji = fs_2, fs_jpim1 
    517                zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    518                zdy2 = e2v(ji,jj) * e2v(ji,jj) 
     569         DO jl = 1, ipl 
     570            DO jj = 1, jpjm1 
     571               DO ji = fs_2, fs_jpim1 
     572                  zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     573                  zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    519574!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj) 
    520                zdy4 = zdy2 * zdy2 
    521                pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                                 ( pt  (ji,jj+1) + pt  (ji,jj)      & 
    522                   &                                                     -          zcv * ( pt  (ji,jj+1) - pt  (ji,jj) ) )  & 
    523                   &        + z1_6   * zdy2 * ( zcv*zcv - 1._wp ) *     (                   ztv2(ji,jj+1) + ztv2(ji,jj)      & 
    524                   &                                                     - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) )  & 
    525                   &        + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1) + ztv4(ji,jj)      & 
    526                   &                                               - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1) - ztv4(ji,jj) ) ) ) 
     575                  zdy4 = zdy2 * zdy2 
     576                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)      & 
     577                     &                                                     -          zcv * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) )  & 
     578                     &        + z1_6   * zdy2 * ( zcv*zcv - 1._wp ) *     (                   ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)      & 
     579                     &                                                     - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) )  & 
     580                     &        + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl)      & 
     581                     &                                               - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) ) 
     582               END DO 
    527583            END DO 
    528584         END DO 
     
    533589    
    534590   
    535    SUBROUTINE nonosc_2d( pbef, paa, pbb, paft, pdt ) 
     591   SUBROUTINE nonosc_2d( ipl, pbef, paa, pbb, paft, pdt ) 
    536592      !!--------------------------------------------------------------------- 
    537593      !!                    ***  ROUTINE nonosc  *** 
     
    546602      !!       in-space based differencing for fluid 
    547603      !!---------------------------------------------------------------------- 
    548       REAL(wp)                     , INTENT(in   ) ::   pdt          ! tracer time-step 
    549       REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pbef, paft   ! before & after field 
    550       REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) ::   paa, pbb     ! monotonic fluxes in the 2 directions 
    551       ! 
    552       INTEGER  ::   ji, jj    ! dummy loop indices 
     604      INTEGER                          , INTENT(in   ) ::   ipl          ! third dimension of tracer array 
     605      REAL(wp)                         , INTENT(in   ) ::   pdt          ! tracer time-step 
     606      REAL(wp), DIMENSION (jpi,jpj,ipl), INTENT(in   ) ::   pbef, paft   ! before & after field 
     607      REAL(wp), DIMENSION (jpi,jpj,ipl), INTENT(inout) ::   paa, pbb     ! monotonic fluxes in the 2 directions 
     608      ! 
     609      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
    553610      INTEGER  ::   ikm1      ! local integer 
    554611      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zsml, z1_dt   ! local scalars 
    555612      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      - 
    556       REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zmsk, zdiv 
     613      REAL(wp), DIMENSION(jpi,jpj) :: zbup, zbdo, zmsk 
     614      REAL(wp), DIMENSION(jpi,jpj,ipl) :: zbetup, zbetdo, zdiv 
    557615      !!---------------------------------------------------------------------- 
    558616      ! 
     
    561619 
    562620      ! test on divergence 
    563       DO jj = 2, jpjm1 
    564          DO ji = fs_2, fs_jpim1   ! vector opt.   
    565             zdiv(ji,jj) =  - (  paa(ji,jj) - paa(ji-1,jj  )   & 
    566                &              + pbb(ji,jj) - pbb(ji  ,jj-1) )   
     621      DO jl = 1, ipl 
     622         DO jj = 2, jpjm1 
     623            DO ji = fs_2, fs_jpim1   ! vector opt.   
     624               zdiv(ji,jj,jl) =  - (  paa(ji,jj,jl) - paa(ji-1,jj  ,jl)   & 
     625                  &              +    pbb(ji,jj,jl) - pbb(ji  ,jj-1,jl) )   
     626            END DO 
    567627         END DO 
    568628      END DO 
    569629      CALL lbc_lnk( 'icedyn_adv_umx', zdiv, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign) 
    570630 
    571       ! Determine ice masks for before and after tracers  
    572       WHERE( pbef(:,:) == 0._wp .AND. paft(:,:) == 0._wp .AND. zdiv(:,:) == 0._wp )   ;   zmsk(:,:) = 0._wp 
    573       ELSEWHERE                                                                       ;   zmsk(:,:) = 1._wp * tmask(:,:,1) 
    574       END WHERE 
    575  
    576       ! Search local extrema 
    577       ! -------------------- 
    578       ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land 
    579 !      zbup(:,:) = MAX( pbef(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) ),   & 
    580 !         &             paft(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) )  ) 
    581 !      zbdo(:,:) = MIN( pbef(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) ),   & 
    582 !         &             paft(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) )  ) 
    583       zbup(:,:) = MAX( pbef(:,:) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) ),   & 
    584          &             paft(:,:) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) )  ) 
    585       zbdo(:,:) = MIN( pbef(:,:) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) ),   & 
    586          &             paft(:,:) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) )  ) 
     631      DO jl = 1, ipl 
     632         ! Determine ice masks for before and after tracers  
     633         WHERE( pbef(:,:,jl) == 0._wp .AND. paft(:,:,jl) == 0._wp .AND. zdiv(:,:,jl) == 0._wp )  
     634            zmsk(:,:) = 0._wp 
     635         ELSEWHERE                                                                                    
     636            zmsk(:,:) = 1._wp * tmask(:,:,1) 
     637         END WHERE 
     638 
     639         ! Search local extrema 
     640         ! -------------------- 
     641         ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land 
     642!         zbup(:,:) = MAX( pbef(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) ),   & 
     643!            &             paft(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) )  ) 
     644!         zbdo(:,:) = MIN( pbef(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) ),   & 
     645!            &             paft(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) )  ) 
     646         zbup(:,:) = MAX( pbef(:,:,jl) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) ),   & 
     647            &             paft(:,:,jl) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) )  ) 
     648         zbdo(:,:) = MIN( pbef(:,:,jl) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) ),   & 
     649            &             paft(:,:,jl) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) )  ) 
    587650 
    588651      z1_dt = 1._wp / pdt 
     
    595658               &                       zbdo(ji  ,jj-1), zbdo(ji  ,jj+1)    ) 
    596659               ! 
    597             zpos = MAX( 0., paa(ji-1,jj  ) ) - MIN( 0., paa(ji  ,jj  ) )   &        ! positive/negative  part of the flux 
    598                & + MAX( 0., pbb(ji  ,jj-1) ) - MIN( 0., pbb(ji  ,jj  ) ) 
    599             zneg = MAX( 0., paa(ji  ,jj  ) ) - MIN( 0., paa(ji-1,jj  ) )   & 
    600                & + MAX( 0., pbb(ji  ,jj  ) ) - MIN( 0., pbb(ji  ,jj-1) ) 
     660               zpos = MAX( 0., paa(ji-1,jj  ,jl) ) - MIN( 0., paa(ji  ,jj  ,jl) )   &        ! positive/negative  part of the flux 
     661                  & + MAX( 0., pbb(ji  ,jj-1,jl) ) - MIN( 0., pbb(ji  ,jj  ,jl) ) 
     662               zneg = MAX( 0., paa(ji  ,jj  ,jl) ) - MIN( 0., paa(ji-1,jj  ,jl) )   & 
     663                  & + MAX( 0., pbb(ji  ,jj  ,jl) ) - MIN( 0., pbb(ji  ,jj-1,jl) ) 
    601664               ! 
    602             zbt = e1e2t(ji,jj) * z1_dt                                   ! up & down beta terms 
    603             zbetup(ji,jj) = ( zup         - paft(ji,jj) ) / ( zpos + zsml ) * zbt 
    604             zbetdo(ji,jj) = ( paft(ji,jj) - zdo         ) / ( zneg + zsml ) * zbt 
     665               zbt = e1e2t(ji,jj) * z1_dt                                   ! up & down beta terms 
     666               zbetup(ji,jj,jl) = ( zup            - paft(ji,jj,jl) ) / ( zpos + zsml ) * zbt 
     667               zbetdo(ji,jj,jl) = ( paft(ji,jj,jl) - zdo         )    / ( zneg + zsml ) * zbt 
     668            END DO 
    605669         END DO 
    606670      END DO 
    607671      CALL lbc_lnk_multi( 'icedyn_adv_umx', zbetup, 'T', 1., zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign) 
    608672 
    609       ! monotonic flux in the i & j direction (paa & pbb) 
    610       ! ------------------------------------- 
    611       DO jj = 2, jpjm1 
    612          DO ji = 1, fs_jpim1   ! vector opt. 
    613             zau = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji+1,jj) ) 
    614             zbu = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji+1,jj) ) 
    615             zcu = 0.5  + SIGN( 0.5 , paa(ji,jj) ) 
    616             ! 
    617             paa(ji,jj) = paa(ji,jj) * ( zcu * zau + ( 1._wp - zcu) * zbu ) 
    618          END DO 
    619       END DO 
    620       ! 
    621       DO jj = 1, jpjm1 
    622          DO ji = fs_2, fs_jpim1   ! vector opt. 
    623             zav = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji,jj+1) ) 
    624             zbv = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji,jj+1) ) 
    625             zcv = 0.5  + SIGN( 0.5 , pbb(ji,jj) ) 
    626             ! 
    627             pbb(ji,jj) = pbb(ji,jj) * ( zcv * zav + ( 1._wp - zcv) * zbv ) 
     673      DO jl = 1, ipl 
     674         ! monotonic flux in the i & j direction (paa & pbb) 
     675         ! ------------------------------------- 
     676         DO jj = 2, jpjm1 
     677            DO ji = 1, fs_jpim1   ! vector opt. 
     678               zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) ) 
     679               zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) ) 
     680               zcu = 0.5  + SIGN( 0.5 , paa(ji,jj,jl) ) 
     681               ! 
     682               paa(ji,jj,jl) = paa(ji,jj,jl) * ( zcu * zau + ( 1._wp - zcu) * zbu ) 
     683            END DO 
     684         END DO 
     685         ! 
     686         DO jj = 1, jpjm1 
     687            DO ji = fs_2, fs_jpim1   ! vector opt. 
     688               zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) ) 
     689               zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) ) 
     690               zcv = 0.5  + SIGN( 0.5 , pbb(ji,jj,jl) ) 
     691               ! 
     692               pbb(ji,jj,jl) = pbb(ji,jj,jl) * ( zcv * zav + ( 1._wp - zcv) * zbv ) 
     693            END DO 
    628694         END DO 
    629695      END DO 
Note: See TracChangeset for help on using the changeset viewer.