Changeset 10180


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

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

Location:
NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src
Files:
4 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 
  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icethd.F90

    r10170 r10180  
    212212         END DO 
    213213 
    214          IF( lk_mpp )         CALL mpp_ini_ice( npti , numout ) 
    215  
    216214         IF( npti > 0 ) THEN  ! If there is no ice, do nothing. 
    217215            !                                                                 
     
    250248            !                                                       ! --- & Move to 2D arrays --- ! 
    251249            ! 
    252             IF( lk_mpp )      CALL mpp_comm_free( ncomm_ice ) !RB necessary ?? 
    253250         ENDIF 
    254251         ! 
  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icethd_zdf_bl99.F90

    r10069 r10180  
    8383      INTEGER, DIMENSION(jpij) ::   jm_min    ! reference number of top equation 
    8484      INTEGER, DIMENSION(jpij) ::   jm_max    ! reference number of bottom equation 
    85       ! 
     85 
     86      LOGICAL, DIMENSION(jpij) ::   l_T_converged   ! true when T converges (per grid point) 
     87! 
    8688      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system 
    8789      REAL(wp) ::   zg1       =  2._wp        ! 
     
    113115      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsb        ! Temporary temperature in the snow to check the convergence 
    114116      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i    ! Ice thermal conductivity 
     117      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i_cp ! copy 
    115118      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradtr_i    ! Radiation transmitted through the ice 
    116119      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradab_i    ! Radiation absorbed in the ice 
     
    201204      ! 
    202205      iconv    = 0          ! number of iterations 
    203       zdti_max = 1000._wp   ! maximal value of error on all points 
    204       ! 
     206      ! 
     207      l_T_converged(:) = .FALSE. 
    205208      !                                                          !============================! 
    206       DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max )   ! Iterative procedure begins ! 
     209      ! Convergence calculated until all sub-domain grid points have converged 
     210      ! Calculations keep going for all grid points until sub-domain convergence (vectorisation optimisation) 
     211      ! but values are not taken into account (results independant of MPI partitioning) 
     212      ! 
     213      DO WHILE ( ( .NOT. ALL (l_T_converged(1:npti)) ) .AND. iconv < iconv_max )   ! Iterative procedure begins ! 
    207214         !                                                       !============================! 
    208215         iconv = iconv + 1 
     
    217224            ! 
    218225            DO ji = 1, npti 
    219                ztcond_i(ji,0)      = rcnd_i + zbeta * sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 
    220                ztcond_i(ji,nlay_i) = rcnd_i + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 ) 
     226               ztcond_i_cp(ji,0)      = rcnd_i + zbeta * sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 
     227               ztcond_i_cp(ji,nlay_i) = rcnd_i + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 ) 
    221228            END DO 
    222229            DO jk = 1, nlay_i-1 
    223230               DO ji = 1, npti 
    224                   ztcond_i(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /  & 
    225                      &                       MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) 
     231                  ztcond_i_cp(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /  & 
     232                     &                         MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) 
    226233               END DO 
    227234            END DO 
     
    230237            ! 
    231238            DO ji = 1, npti 
    232                ztcond_i(ji,0)      = rcnd_i + 0.09_wp  *  sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  & 
    233                   &                         - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 
    234                ztcond_i(ji,nlay_i) = rcnd_i + 0.09_wp  *  sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )  & 
    235                   &                         - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 
     239               ztcond_i_cp(ji,0)      = rcnd_i + 0.09_wp  *  sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  & 
     240                  &                           - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 
     241               ztcond_i_cp(ji,nlay_i) = rcnd_i + 0.09_wp  *  sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )  & 
     242                  &                           - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 
    236243            END DO 
    237244            DO jk = 1, nlay_i-1 
    238245               DO ji = 1, npti 
    239                   ztcond_i(ji,jk) = rcnd_i + 0.09_wp  *   0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /        & 
    240                      &                      MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )  & 
    241                      &                     - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) 
     246                  ztcond_i_cp(ji,jk) = rcnd_i + 0.09_wp  *   0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /        & 
     247                     &                        MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )  & 
     248                     &                       - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) 
    242249               END DO 
    243250            END DO 
    244251            ! 
    245252         ENDIF 
    246          ztcond_i(1:npti,:) = MAX( zkimin, ztcond_i(1:npti,:) )         
     253 
     254         ! Variable used after iterations 
     255         ! Value must be frozen after convergence for MPP independance reason 
     256         DO ji = 1, npti 
     257            IF ( .NOT. l_T_converged(ji) ) & 
     258               ztcond_i(ji,:) = MAX( zkimin, ztcond_i_cp(ji,:) )         
     259         END DO 
    247260         ! 
    248261         !--- G(he) : enhancement of thermal conductivity in mono-category case 
     
    270283         !----------------- 
    271284         !--- Snow 
     285         ! Variable used after iterations 
     286         ! Value must be frozen after convergence for MPP independance reason 
    272287         DO jk = 0, nlay_s-1 
    273288            DO ji = 1, npti 
    274                zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji) 
     289               IF ( .NOT. l_T_converged(ji) ) & 
     290                  zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji) 
    275291            END DO 
    276292         END DO 
    277293         DO ji = 1, npti   ! Snow-ice interface 
    278             zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) 
    279             IF( zfac > epsi10 ) THEN 
    280                zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac 
    281             ELSE 
    282                zkappa_s(ji,nlay_s) = 0._wp 
     294            IF ( .NOT. l_T_converged(ji) ) THEN 
     295               zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) 
     296               IF( zfac > epsi10 ) THEN 
     297                  zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac 
     298               ELSE 
     299                  zkappa_s(ji,nlay_s) = 0._wp 
     300               ENDIF 
    283301            ENDIF 
    284302         END DO 
    285303 
    286304         !--- Ice 
     305         ! Variable used after iterations 
     306         ! Value must be frozen after convergence for MPP independance reason 
    287307         DO jk = 0, nlay_i 
    288308            DO ji = 1, npti 
    289                zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji) 
     309               IF ( .NOT. l_T_converged(ji) ) & 
     310                  zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji) 
    290311            END DO 
    291312         END DO 
    292313         DO ji = 1, npti   ! Snow-ice interface 
    293             zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 
     314            IF ( .NOT. l_T_converged(ji) ) & 
     315               zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 
    294316         END DO 
    295317         ! 
     
    326348            ! update of the non solar flux according to the update in T_su 
    327349            DO ji = 1, npti 
    328                qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 
     350               ! Variable used after iterations 
     351               ! Value must be frozen after convergence for MPP independance reason 
     352               IF ( .NOT. l_T_converged(ji) ) & 
     353                  qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 
    329354            END DO 
    330355 
     
    496521            ! ice temperatures 
    497522            DO ji = 1, npti 
    498                t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 
     523               ! Variable used after iterations 
     524               ! Value must be frozen after convergence for MPP independance reason 
     525               IF ( .NOT. l_T_converged(ji) ) & 
     526                  t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 
    499527            END DO 
    500528 
     
    502530               DO ji = 1, npti 
    503531                  jk = jm - nlay_s - 1 
    504                   t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
     532                  IF ( .NOT. l_T_converged(ji) ) & 
     533                     t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
    505534               END DO 
    506535            END DO 
    507536 
    508537            DO ji = 1, npti 
    509                ! snow temperatures       
    510                IF( h_s_1d(ji) > 0._wp ) THEN 
    511                   t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 
    512                ENDIF 
    513                ! surface temperature 
    514                ztsub(ji) = t_su_1d(ji) 
    515                IF( t_su_1d(ji) < rt0 ) THEN 
    516                   t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  & 
    517                      &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 
     538               ! Variables used after iterations 
     539               ! Value must be frozen after convergence for MPP independance reason 
     540               IF ( .NOT. l_T_converged(ji) ) THEN 
     541                  ! snow temperatures       
     542                  IF( h_s_1d(ji) > 0._wp ) THEN 
     543                     t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 
     544                  ENDIF 
     545                  ! surface temperature 
     546                  ztsub(ji) = t_su_1d(ji) 
     547                  IF( t_su_1d(ji) < rt0 ) THEN 
     548                     t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  & 
     549                        &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 
     550                  ENDIF 
    518551               ENDIF 
    519552            END DO 
     
    524557            ! check that nowhere it has started to melt 
    525558            ! zdti_max is a measure of error, it has to be under zdti_bnd 
    526             zdti_max = 0._wp 
    527             DO ji = 1, npti 
    528                t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 
    529                zdti_max    = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )      
    530             END DO 
    531  
    532             DO jk = 1, nlay_s 
    533                DO ji = 1, npti 
    534                   t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 
    535                   zdti_max      = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 
    536                END DO 
    537             END DO 
    538  
    539             DO jk = 1, nlay_i 
    540                DO ji = 1, npti 
    541                   ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0  
    542                   t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 
    543                   zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 
    544                END DO 
    545             END DO 
    546  
    547             ! Compute spatial maximum over all errors 
    548             ! note that this could be optimized substantially by iterating only the non-converging points 
    549             IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 
    550          ! 
     559 
     560            DO ji = 1, npti 
     561 
     562               zdti_max = 0._wp 
     563 
     564               IF ( .NOT. l_T_converged(ji) ) THEN 
     565                  t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 
     566                  zdti_max    = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) ) 
     567 
     568                  t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp ) 
     569                  zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) ) 
     570 
     571                  DO jk = 1, nlay_i 
     572                     ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0 
     573                     t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 
     574                     zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 
     575                  END DO 
     576 
     577                  IF ( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE. 
     578 
     579               ENDIF 
     580 
     581            END DO 
     582 
    551583         !----------------------------------------! 
    552584         !                                        ! 
     
    670702             
    671703            ! ice temperatures 
    672            DO ji = 1, npti 
    673                t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 
     704            DO ji = 1, npti 
     705               ! Variable used after iterations 
     706               ! Value must be frozen after convergence for MPP independance reason 
     707               IF ( .NOT. l_T_converged(ji) ) & 
     708                  t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 
    674709            END DO 
    675710 
    676711            DO jm = nlay_i + nlay_s, nlay_s + 2, -1 
    677712               DO ji = 1, npti 
    678                   jk = jm - nlay_s - 1 
    679                   t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
     713                  IF ( .NOT. l_T_converged(ji) ) THEN 
     714                     jk = jm - nlay_s - 1 
     715                     t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
     716                  ENDIF 
    680717               END DO 
    681718            END DO 
     
    683720            ! snow temperatures       
    684721            DO ji = 1, npti 
    685                IF( h_s_1d(ji) > 0._wp ) THEN 
    686                   t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 
     722               ! Variable used after iterations 
     723               ! Value must be frozen after convergence for MPP independance reason 
     724               IF ( .NOT. l_T_converged(ji) ) THEN 
     725                  IF( h_s_1d(ji) > 0._wp ) THEN 
     726                     t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 
     727                  ENDIF 
    687728               ENDIF 
    688729            END DO 
     
    693734            ! check that nowhere it has started to melt 
    694735            ! zdti_max is a measure of error, it has to be under zdti_bnd 
    695             zdti_max = 0._wp 
    696  
    697             DO jk = 1, nlay_s 
    698                DO ji = 1, npti 
    699                   t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 
    700                   zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 
    701                END DO 
    702             END DO 
    703              
    704             DO jk = 1, nlay_i 
    705                DO ji = 1, npti 
    706                   ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0  
    707                   t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 
    708                   zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 
    709                END DO 
    710             END DO 
    711  
    712             ! Compute spatial maximum over all errors 
    713             ! note that this could be optimized substantially by iterating only the non-converging points 
    714             IF( lk_mpp )   CALL mpp_max( zdti_max, kcom=ncomm_ice ) 
     736 
     737            DO ji = 1, npti 
     738 
     739               zdti_max = 0._wp 
     740 
     741               IF ( .NOT. l_T_converged(ji) ) THEN 
     742                  ! t_s 
     743                  t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp ) 
     744                  zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) ) 
     745                  ! t_i 
     746                  DO jk = 0, nlay_i 
     747                     ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0  
     748                     t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 
     749                     zdti_max      =  MAX ( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 
     750                  END DO 
     751 
     752                  IF ( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE. 
     753 
     754               ENDIF 
     755 
     756            END DO 
    715757 
    716758         ENDIF ! k_jules 
  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/LBC/lib_mpp.F90

    r10172 r10180  
    8585   PUBLIC   mpp_max_multiple 
    8686   PUBLIC   mppscatter, mppgather 
    87    PUBLIC   mpp_ini_ice, mpp_ini_znl 
     87   PUBLIC   mpp_ini_znl 
    8888   PUBLIC   mppsize 
    8989   PUBLIC   mppsend, mpprecv                          ! needed by TAM and ICB routines 
     
    133133 
    134134   INTEGER :: MPI_SUMDD 
    135  
    136    ! variables used in case of sea-ice 
    137    INTEGER, PUBLIC ::   ncomm_ice       !: communicator made by the processors with sea-ice (public so that it can be freed in icethd) 
    138    INTEGER         ::   ngrp_iworld     !  group ID for the world processors (for rheology) 
    139    INTEGER         ::   ngrp_ice        !  group ID for the ice processors (for rheology) 
    140    INTEGER         ::   ndim_rank_ice   !  number of 'ice' processors 
    141    INTEGER         ::   n_ice_root      !  number (in the comm_ice) of proc 0 in the ice comm 
    142    INTEGER, DIMENSION(:), ALLOCATABLE, SAVE ::   nrank_ice     ! dimension ndim_rank_ice 
    143135 
    144136   ! variables used for zonal integration 
     
    10111003 
    10121004 
    1013    SUBROUTINE mpp_ini_ice( pindic, kumout ) 
    1014       !!---------------------------------------------------------------------- 
    1015       !!               ***  routine mpp_ini_ice  *** 
    1016       !! 
    1017       !! ** Purpose :   Initialize special communicator for ice areas 
    1018       !!      condition together with global variables needed in the ddmpp folding 
    1019       !! 
    1020       !! ** Method  : - Look for ice processors in ice routines 
    1021       !!              - Put their number in nrank_ice 
    1022       !!              - Create groups for the world processors and the ice processors 
    1023       !!              - Create a communicator for ice processors 
    1024       !! 
    1025       !! ** output 
    1026       !!      njmppmax = njmpp for northern procs 
    1027       !!      ndim_rank_ice = number of processors with ice 
    1028       !!      nrank_ice (ndim_rank_ice) = ice processors 
    1029       !!      ngrp_iworld = group ID for the world processors 
    1030       !!      ngrp_ice = group ID for the ice processors 
    1031       !!      ncomm_ice = communicator for the ice procs. 
    1032       !!      n_ice_root = number (in the world) of proc 0 in the ice comm. 
    1033       !! 
    1034       !!---------------------------------------------------------------------- 
    1035       INTEGER, INTENT(in) ::   pindic 
    1036       INTEGER, INTENT(in) ::   kumout   ! ocean.output logical unit 
    1037       !! 
    1038       INTEGER :: jjproc 
    1039       INTEGER :: ii, ierr 
    1040       INTEGER, ALLOCATABLE, DIMENSION(:) ::   kice 
    1041       INTEGER, ALLOCATABLE, DIMENSION(:) ::   zwork 
    1042       !!---------------------------------------------------------------------- 
    1043       ! 
    1044       ALLOCATE( kice(jpnij), zwork(jpnij), STAT=ierr ) 
    1045       IF( ierr /= 0 ) THEN 
    1046          WRITE(kumout, cform_err) 
    1047          WRITE(kumout,*) 'mpp_ini_ice : failed to allocate 2, 1D arrays (jpnij in length)' 
    1048          CALL mppstop 
    1049       ENDIF 
    1050  
    1051       ! Look for how many procs with sea-ice 
    1052       ! 
    1053       kice = 0 
    1054       DO jjproc = 1, jpnij 
    1055          IF( jjproc == narea .AND. pindic .GT. 0 )   kice(jjproc) = 1 
    1056       END DO 
    1057       ! 
    1058       zwork = 0 
    1059       CALL MPI_ALLREDUCE( kice, zwork, jpnij, mpi_integer, mpi_sum, mpi_comm_oce, ierr ) 
    1060       ndim_rank_ice = SUM( zwork ) 
    1061  
    1062       ! Allocate the right size to nrank_north 
    1063       IF( ALLOCATED ( nrank_ice ) )   DEALLOCATE( nrank_ice ) 
    1064       ALLOCATE( nrank_ice(ndim_rank_ice) ) 
    1065       ! 
    1066       ii = 0 
    1067       nrank_ice = 0 
    1068       DO jjproc = 1, jpnij 
    1069          IF( zwork(jjproc) == 1) THEN 
    1070             ii = ii + 1 
    1071             nrank_ice(ii) = jjproc -1 
    1072          ENDIF 
    1073       END DO 
    1074  
    1075       ! Create the world group 
    1076       CALL MPI_COMM_GROUP( mpi_comm_oce, ngrp_iworld, ierr ) 
    1077  
    1078       ! Create the ice group from the world group 
    1079       CALL MPI_GROUP_INCL( ngrp_iworld, ndim_rank_ice, nrank_ice, ngrp_ice, ierr ) 
    1080  
    1081       ! Create the ice communicator , ie the pool of procs with sea-ice 
    1082       CALL MPI_COMM_CREATE( mpi_comm_oce, ngrp_ice, ncomm_ice, ierr ) 
    1083  
    1084       ! Find proc number in the world of proc 0 in the north 
    1085       ! The following line seems to be useless, we just comment & keep it as reminder 
    1086       ! CALL MPI_GROUP_TRANSLATE_RANKS(ngrp_ice,1,0,ngrp_iworld,n_ice_root,ierr) 
    1087       ! 
    1088       CALL MPI_GROUP_FREE(ngrp_ice, ierr) 
    1089       CALL MPI_GROUP_FREE(ngrp_iworld, ierr) 
    1090  
    1091       DEALLOCATE(kice, zwork) 
    1092       ! 
    1093    END SUBROUTINE mpp_ini_ice 
    1094  
    1095  
    10961005   SUBROUTINE mpp_ini_znl( kumout ) 
    10971006      !!---------------------------------------------------------------------- 
     
    16621571   LOGICAL, PUBLIC, PARAMETER ::   lk_mpp = .FALSE.      !: mpp flag 
    16631572   LOGICAL, PUBLIC            ::   ln_nnogather          !: namelist control of northfold comms (needed here in case "key_mpp_mpi" is not used) 
    1664    INTEGER :: ncomm_ice 
    16651573   INTEGER, PUBLIC            ::   mpi_comm_oce          ! opa local communicator 
    16661574   !!---------------------------------------------------------------------- 
     
    18151723      STOP      ! non MPP case, just stop the run 
    18161724   END SUBROUTINE mppstop 
    1817  
    1818    SUBROUTINE mpp_ini_ice( kcom, knum ) 
    1819       INTEGER :: kcom, knum 
    1820       WRITE(*,*) 'mpp_ini_ice: You should not have seen this print! error?', kcom, knum 
    1821    END SUBROUTINE mpp_ini_ice 
    18221725 
    18231726   SUBROUTINE mpp_ini_znl( knum ) 
Note: See TracChangeset for help on using the changeset viewer.