Changeset 11612


Ignore:
Timestamp:
2019-09-29T10:55:44+02:00 (13 months ago)
Author:
clem
Message:

make Prather the default advection scheme for ice simulations and drastically reduce the number of communications in Prather routine (by a factor of 10)

Location:
NEMO/trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/cfgs/SHARED/namelist_ice_ref

    r11601 r11612  
    103103&namdyn_adv     !   Ice advection 
    104104!------------------------------------------------------------------------------ 
    105    ln_adv_Pra       = .false.         !  Advection scheme (Prather) 
    106    ln_adv_UMx       = .true.          !  Advection scheme (Ultimate-Macho) 
     105   ln_adv_Pra       = .true.         !  Advection scheme (Prather) 
     106   ln_adv_UMx       = .false.          !  Advection scheme (Ultimate-Macho) 
    107107      nn_UMx        =   5             !     order of the scheme for UMx (1-5 ; 20=centered 2nd order) 
    108108/ 
  • NEMO/trunk/cfgs/SPITZ12/EXPREF/namelist_ice_cfg

    r11548 r11612  
    5050&namdyn_adv     !   Ice advection 
    5151!------------------------------------------------------------------------------ 
     52   ln_adv_Pra       = .false.         !  Advection scheme (Prather) 
     53   ln_adv_UMx       = .true.          !  Advection scheme (Ultimate-Macho) 
     54      nn_UMx        =   5             !     order of the scheme for UMx (1-5 ; 20=centered 2nd order) 
    5255/ 
    5356!------------------------------------------------------------------------------ 
  • NEMO/trunk/src/ICE/icectl.F90

    r11601 r11612  
    158158            !    only check for Prather because Ultimate-Macho uses corrective fluxes (wfx etc) 
    159159            !    so the formulation for conservation is different (and not coded)  
    160             !    it does not mean UM is not conservative (it is checked with above prints) 
    161             IF( ln_adv_Pra .AND. ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 
    162                &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice 
     160            !    it does not mean UM is not conservative (it is checked with above prints) => update (09/2019): same for Prather now 
     161            !IF( ln_adv_Pra .AND. ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 
     162            !   &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice 
    163163         ENDIF 
    164164         ! 
  • NEMO/trunk/src/ICE/icedyn_adv_pra.F90

    r10425 r11612  
    1919   USE ice            ! sea-ice variables 
    2020   USE sbc_oce , ONLY : nn_fsbc   ! frequency of sea-ice call 
     21   USE icevar         ! sea-ice: operations 
    2122   ! 
    2223   USE in_out_manager ! I/O manager 
     
    2526   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero) 
    2627   USE lbclnk         ! lateral boundary conditions (or mpp links) 
    27    USE prtctl         ! Print control 
    2828 
    2929   IMPLICIT NONE 
     
    3939   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxsal, sysal, sxxsal, syysal, sxysal   ! ice salinity 
    4040   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxage, syage, sxxage, syyage, sxyage   ! ice age 
    41    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:   ::   sxopw, syopw, sxxopw, syyopw, sxyopw   ! open water in sea ice 
     41   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxopw, syopw, sxxopw, syyopw, sxyopw   ! open water in sea ice 
    4242   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sxc0 , syc0 , sxxc0 , syyc0 , sxyc0    ! snow layers heat content 
    4343   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sxe  , sye  , sxxe  , syye  , sxye     ! ice layers heat content 
     
    8282      ! 
    8383      INTEGER  ::   jk, jl, jt              ! dummy loop indices 
    84       INTEGER  ::   initad                  ! number of sub-timestep for the advection 
    85       REAL(wp) ::   zcfl , zusnit           !   -      - 
    86       REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   zarea 
    87       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0opw 
    88       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0ice, z0snw, z0ai, z0smi, z0oi 
    89       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0ap , z0vp 
    90       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   z0es 
    91       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   z0ei 
     84      INTEGER  ::   icycle                  ! number of sub-timestep for the advection 
     85      REAL(wp) ::   zdt                     !   -      - 
     86      REAL(wp), DIMENSION(1)                  ::   zcflprv, zcflnow   ! for global communication 
     87      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zarea 
     88      REAL(wp), DIMENSION(jpi,jpj,1)          ::   z0opw 
     89      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ice, z0snw, z0ai, z0smi, z0oi 
     90      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp 
     91      REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   z0es 
     92      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   z0ei 
    9293      !!---------------------------------------------------------------------- 
    9394      ! 
    9495      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme' 
    9596      ! 
    96       ALLOCATE( zarea(jpi,jpj)    , z0opw(jpi,jpj, 1 ) , z0ice(jpi,jpj,jpl) , z0snw(jpi,jpj,jpl) ,                       & 
    97          &      z0ai(jpi,jpj,jpl) , z0smi(jpi,jpj,jpl) , z0oi (jpi,jpj,jpl) , z0ap (jpi,jpj,jpl) , z0vp(jpi,jpj,jpl) ,   & 
    98          &      z0es (jpi,jpj,nlay_s,jpl), z0ei(jpi,jpj,nlay_i,jpl) ) 
    99       ! 
    100       ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !         
    101       zcfl  =            MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) 
    102       zcfl  = MAX( zcfl, MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 
    103       CALL mpp_max( 'icedyn_adv_pra', zcfl ) 
     97      ! --- If ice drift is too fast, use  subtime steps for advection (CFL test for stability) --- ! 
     98      !        Note: the advection split is applied at the next time-step in order to avoid blocking global comm. 
     99      !              this should not affect too much the stability 
     100      zcflnow(1) =                  MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) 
     101      zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 
    104102       
    105       IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp 
    106       ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp 
     103      ! non-blocking global communication send zcflnow and receive zcflprv 
     104      CALL mpp_delay_max( 'icedyn_adv_pra', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 ) 
     105 
     106      IF( zcflprv(1) > .5 ) THEN   ;   icycle = 2 
     107      ELSE                         ;   icycle = 1 
    107108      ENDIF 
     109      zdt = rdt_ice / REAL(icycle) 
    108110       
    109       zarea(:,:) = e1e2t(:,:) 
    110111      !------------------------- 
    111112      ! transported fields                                         
     
    113114      z0opw(:,:,1) = pato_i(:,:) * e1e2t(:,:)              ! Open water area  
    114115      DO jl = 1, jpl 
    115          z0snw(:,:,jl) = pv_s (:,:,  jl) * e1e2t(:,:)     ! Snow volume 
    116          z0ice(:,:,jl) = pv_i (:,:,  jl) * e1e2t(:,:)     ! Ice  volume 
    117          z0ai (:,:,jl) = pa_i (:,:,  jl) * e1e2t(:,:)     ! Ice area 
    118          z0smi(:,:,jl) = psv_i(:,:,  jl) * e1e2t(:,:)     ! Salt content 
    119          z0oi (:,:,jl) = poa_i(:,:,  jl) * e1e2t(:,:)     ! Age content 
     116         zarea(:,:,jl) = e1e2t(:,:) 
     117         z0snw(:,:,jl) = pv_s (:,:,jl) * e1e2t(:,:)        ! Snow volume 
     118         z0ice(:,:,jl) = pv_i (:,:,jl) * e1e2t(:,:)        ! Ice  volume 
     119         z0ai (:,:,jl) = pa_i (:,:,jl) * e1e2t(:,:)        ! Ice area 
     120         z0smi(:,:,jl) = psv_i(:,:,jl) * e1e2t(:,:)        ! Salt content 
     121         z0oi (:,:,jl) = poa_i(:,:,jl) * e1e2t(:,:)        ! Age content 
    120122         DO jk = 1, nlay_s 
    121123            z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content 
     
    133135      IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==! 
    134136         !                                                 !--------------------------------------------! 
    135          DO jt = 1, initad 
    136             CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area 
    137                &                                      sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    138             CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   & 
    139                &                                      sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    140             DO jl = 1, jpl 
    141                CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
    142                   &                                      sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    143                CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   & 
    144                   &                                      sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    145                CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
    146                   &                                      sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    147                CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   & 
    148                   &                                      sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    149                CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
    150                   &                                      sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    151                CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   & 
    152                   &                                      sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    153                CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---      
    154                   &                                      sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    155                CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   & 
    156                   &                                      sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    157                CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations --- 
    158                   &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    159                CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &  
    160                   &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    161                DO jk = 1, nlay_s                                                               !--- snow heat contents --- 
    162                   CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl),   & 
    163                      &                                      sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) ) 
    164                   CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl),   & 
    165                      &                                      sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) ) 
    166                END DO 
    167                DO jk = 1, nlay_i                                                               !--- ice heat contents --- 
    168                   CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl),   &  
    169                      &                                      sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    170                   CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl),   &  
    171                      &                                      sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    172                END DO 
    173                IF ( ln_pnd_H12 ) THEN 
    174                   CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &    !--- melt pond fraction -- 
    175                      &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  ) 
    176                   CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &  
    177                      &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  ) 
    178                   CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &    !--- melt pond volume   -- 
    179                      &                                      sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  ) 
    180                   CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &  
    181                      &                                      sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  ) 
    182                ENDIF 
    183             END DO 
     137         DO jt = 1, icycle 
     138            CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0opw , sxopw , sxxopw , syopw , syyopw , sxyopw ) !--- open water 
     139            CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0opw , sxopw , sxxopw , syopw , syyopw , sxyopw ) 
     140            CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume 
     141            CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) 
     142            CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) !--- snow volume 
     143            CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) 
     144            CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity 
     145            CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) 
     146            CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) !--- ice concentration 
     147            CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) 
     148            CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) !--- ice age 
     149            CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) 
     150            ! 
     151            DO jk = 1, nlay_s                                                                             !--- snow heat content 
     152               CALL adv_x( zdt, pu_ice, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   & 
     153                  &                                   sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 
     154               CALL adv_y( zdt, pv_ice, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   & 
     155                  &                                   sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 
     156            END DO 
     157            DO jk = 1, nlay_i                                                                             !--- ice heat content 
     158               CALL adv_x( zdt, pu_ice, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   &  
     159                  &                                   sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
     160               CALL adv_y( zdt, pv_ice, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   &  
     161                  &                                   sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
     162            END DO 
     163            ! 
     164            IF ( ln_pnd_H12 ) THEN 
     165               CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction 
     166               CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )  
     167               CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume 
     168               CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )  
     169            ENDIF 
    184170         END DO 
    185171      !                                                    !--------------------------------------------! 
    186172      ELSE                                                 !== even ice time step:  adv_y then adv_x  ==! 
    187173         !                                                 !--------------------------------------------! 
    188          DO jt = 1, initad 
    189             CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area 
    190                &                                      sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    191             CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   & 
    192                &                                      sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    193             DO jl = 1, jpl 
    194                CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
    195                   &                                      sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    196                CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   & 
    197                   &                                      sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    198                CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
    199                   &                                      sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    200                CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   & 
    201                   &                                      sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    202                CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
    203                   &                                      sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    204                CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   & 
    205                   &                                      sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    206                CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &   !--- ice age      --- 
    207                   &                                      sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    208                CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   & 
    209                   &                                      sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    210                CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations --- 
    211                   &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    212                CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
    213                   &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    214                DO jk = 1, nlay_s                                                             !--- snow heat contents --- 
    215                   CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl),   & 
    216                      &                                      sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) ) 
    217                   CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl),   & 
    218                      &                                      sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) ) 
    219                END DO 
    220                DO jk = 1, nlay_i                                                             !--- ice heat contents --- 
    221                   CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl),   &  
    222                      &                                      sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    223                   CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl),   &  
    224                      &                                      sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    225                END DO 
    226                IF ( ln_pnd_H12 ) THEN 
    227                   CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &   !--- melt pond fraction --- 
    228                      &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  ) 
    229                   CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   & 
    230                      &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  ) 
    231                   CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &   !--- melt pond volume   --- 
    232                      &                                      sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  ) 
    233                   CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   & 
    234                      &                                      sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  ) 
    235                ENDIF 
    236             END DO 
     174         DO jt = 1, icycle 
     175            CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0opw , sxopw , sxxopw , syopw , syyopw , sxyopw ) !--- open water 
     176            CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0opw , sxopw , sxxopw , syopw , syyopw , sxyopw ) 
     177            CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume 
     178            CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) 
     179            CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) !--- snow volume 
     180            CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) 
     181            CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity 
     182            CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) 
     183            CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) !--- ice concentration 
     184            CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) 
     185            CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) !--- ice age 
     186            CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) 
     187            DO jk = 1, nlay_s                                                                             !--- snow heat content 
     188               CALL adv_y( zdt, pv_ice, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   & 
     189                  &                                   sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 
     190               CALL adv_x( zdt, pu_ice, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   & 
     191                  &                                   sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 
     192            END DO 
     193            DO jk = 1, nlay_i                                                                             !--- ice heat content 
     194               CALL adv_y( zdt, pv_ice, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   &  
     195                  &                                   sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
     196               CALL adv_x( zdt, pu_ice, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   &  
     197                  &                                   sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
     198            END DO 
     199            IF ( ln_pnd_H12 ) THEN 
     200               CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction 
     201               CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 
     202               CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume 
     203               CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 
     204            ENDIF 
    237205         END DO 
    238206      ENDIF 
     
    243211      pato_i(:,:) = z0opw(:,:,1) * r1_e1e2t(:,:) * tmask(:,:,1) 
    244212      DO jl = 1, jpl 
    245          pv_i (:,:,  jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    246          pv_s (:,:,  jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    247          psv_i(:,:,  jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    248          poa_i(:,:,  jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    249          pa_i (:,:,  jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     213         pv_i (:,:,jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     214         pv_s (:,:,jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     215         psv_i(:,:,jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     216         poa_i(:,:,jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     217         pa_i (:,:,jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    250218         DO jk = 1, nlay_s 
    251219            pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     
    255223         END DO 
    256224         IF ( ln_pnd_H12 ) THEN 
    257             pa_ip  (:,:,jl) = z0ap (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    258             pv_ip  (:,:,jl) = z0vp (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     225            pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     226            pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    259227         ENDIF 
    260228      END DO 
    261229      ! 
    262       DEALLOCATE( zarea , z0opw , z0ice, z0snw , z0ai , z0smi , z0oi , z0ap , z0vp , z0es, z0ei ) 
     230      ! --- Ensure non-negative fields --- ! 
     231      ! Remove negative values (conservation is ensured) 
     232      !    (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 
     233      CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
    263234      ! 
    264235      IF( lrst_ice )   CALL adv_pra_rst( 'WRITE', kt )   !* write Prather fields in the restart file 
     
    267238    
    268239    
    269    SUBROUTINE adv_x( pdf, put , pcrh, psm , ps0 ,   & 
     240   SUBROUTINE adv_x( pdt, put , pcrh, psm , ps0 ,   & 
    270241      &              psx, psxx, psy , psyy, psxy ) 
    271242      !!---------------------------------------------------------------------- 
     
    275246      !!                variable on x axis 
    276247      !!---------------------------------------------------------------------- 
    277       REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step 
    278       REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0) 
    279       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s] 
    280       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area 
    281       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected 
    282       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments  
    283       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments 
     248      REAL(wp)                  , INTENT(in   ) ::   pdt                ! the time step 
     249      REAL(wp)                  , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0) 
     250      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s] 
     251      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psm                ! area 
     252      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ps0                ! field to be advected 
     253      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psx , psy          ! 1st moments  
     254      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments 
    284255      !!  
    285       INTEGER  ::   ji, jj                               ! dummy loop indices 
    286       REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp         ! local scalars 
     256      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices 
     257      REAL(wp) ::   zs1max, zslpmax, ztemp               ! local scalars 
    287258      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      - 
    288259      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      - 
     
    291262      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q         !  -      - 
    292263      !----------------------------------------------------------------------- 
    293  
    294       ! Limitation of moments.                                            
    295  
    296       zrdt = rdt_ice * pdf      ! If ice drift field is too fast, use an appropriate time step for advection. 
    297  
    298       DO jj = 1, jpj 
    299          DO ji = 1, jpi 
    300             zslpmax = MAX( 0._wp, ps0(ji,jj) ) 
    301             zs1max  = 1.5 * zslpmax 
    302             zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj) ) ) 
    303             zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      & 
    304                &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) )  ) 
    305             rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
    306  
    307             ps0 (ji,jj) = zslpmax   
    308             psx (ji,jj) = zs1new      * rswitch 
    309             psxx(ji,jj) = zs2new      * rswitch 
    310             psy (ji,jj) = psy (ji,jj) * rswitch 
    311             psyy(ji,jj) = psyy(ji,jj) * rswitch 
    312             psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch 
    313          END DO 
     264      ! 
     265      jcat = SIZE( ps0 , 3 )   ! size of input arrays 
     266      ! 
     267      DO jl = 1, jcat   ! loop on categories 
     268         ! 
     269         ! Limitation of moments.                                            
     270         DO jj = 2, jpjm1 
     271            DO ji = 1, jpi 
     272               !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                      
     273               psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 
     274               ! 
     275               zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 
     276               zs1max  = 1.5 * zslpmax 
     277               zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) ) 
     278               zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      & 
     279                  &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) )  ) 
     280               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
     281 
     282               ps0 (ji,jj,jl) = zslpmax   
     283               psx (ji,jj,jl) = zs1new         * rswitch 
     284               psxx(ji,jj,jl) = zs2new         * rswitch 
     285               psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch 
     286               psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch 
     287               psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 
     288            END DO 
     289         END DO 
     290 
     291         !  Calculate fluxes and moments between boxes i<-->i+1               
     292         DO jj = 2, jpjm1                      !  Flux from i to i+1 WHEN u GT 0  
     293            DO ji = 1, jpi 
     294               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 
     295               zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt * e2u(ji,jj) / psm(ji,jj,jl) 
     296               zalfq        =  zalf * zalf 
     297               zalf1        =  1.0 - zalf 
     298               zalf1q       =  zalf1 * zalf1 
     299               ! 
     300               zfm (ji,jj)  =  zalf  *   psm (ji,jj,jl) 
     301               zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) ) 
     302               zfx (ji,jj)  =  zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) ) 
     303               zfxx(ji,jj)  =  zalf  *   psxx(ji,jj,jl) * zalfq 
     304               zfy (ji,jj)  =  zalf  * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 
     305               zfxy(ji,jj)  =  zalfq *   psxy(ji,jj,jl) 
     306               zfyy(ji,jj)  =  zalf  *   psyy(ji,jj,jl) 
     307 
     308               !  Readjust moments remaining in the box. 
     309               psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj) 
     310               ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj) 
     311               psx (ji,jj,jl)  =  zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) ) 
     312               psxx(ji,jj,jl)  =  zalf1  * zalf1q * psxx(ji,jj,jl) 
     313               psy (ji,jj,jl)  =  psy (ji,jj,jl) - zfy(ji,jj) 
     314               psyy(ji,jj,jl)  =  psyy(ji,jj,jl) - zfyy(ji,jj) 
     315               psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl) 
     316            END DO 
     317         END DO 
     318 
     319         DO jj = 2, jpjm1                      !  Flux from i+1 to i when u LT 0. 
     320            DO ji = 1, fs_jpim1 
     321               zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt * e2u(ji,jj) / psm(ji+1,jj,jl)  
     322               zalg  (ji,jj) = zalf 
     323               zalfq         = zalf * zalf 
     324               zalf1         = 1.0 - zalf 
     325               zalg1 (ji,jj) = zalf1 
     326               zalf1q        = zalf1 * zalf1 
     327               zalg1q(ji,jj) = zalf1q 
     328               ! 
     329               zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji+1,jj,jl) 
     330               zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji+1,jj,jl) & 
     331                  &                                   - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 
     332               zfx   (ji,jj) = zfx (ji,jj) + zalfq * (  psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 
     333               zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji+1,jj,jl) * zalfq 
     334               zfy   (ji,jj) = zfy (ji,jj) + zalf  * (  psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 
     335               zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji+1,jj,jl) 
     336               zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji+1,jj,jl) 
     337            END DO 
     338         END DO 
     339 
     340         DO jj = 2, jpjm1                     !  Readjust moments remaining in the box.  
     341            DO ji = fs_2, fs_jpim1 
     342               zbt  =       zbet(ji-1,jj) 
     343               zbt1 = 1.0 - zbet(ji-1,jj) 
     344               ! 
     345               psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) ) 
     346               ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) ) 
     347               psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) ) 
     348               psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl) 
     349               psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) ) 
     350               psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) ) 
     351               psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl) 
     352            END DO 
     353         END DO 
     354 
     355         !   Put the temporary moments into appropriate neighboring boxes.     
     356         DO jj = 2, jpjm1                     !   Flux from i to i+1 IF u GT 0. 
     357            DO ji = fs_2, fs_jpim1 
     358               zbt  =       zbet(ji-1,jj) 
     359               zbt1 = 1.0 - zbet(ji-1,jj) 
     360               psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl) 
     361               zalf          = zbt * zfm(ji-1,jj) / psm(ji,jj,jl) 
     362               zalf1         = 1.0 - zalf 
     363               ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj) 
     364               ! 
     365               ps0 (ji,jj,jl) =  zbt  * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl) 
     366               psx (ji,jj,jl) =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl) 
     367               psxx(ji,jj,jl) =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl)                             & 
     368                  &                     + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  ) & 
     369                  &            + zbt1 * psxx(ji,jj,jl) 
     370               psxy(ji,jj,jl) =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl)             & 
     371                  &                     + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj,jl) ) )   & 
     372                  &            + zbt1 * psxy(ji,jj,jl) 
     373               psy (ji,jj,jl) =  zbt  * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl) 
     374               psyy(ji,jj,jl) =  zbt  * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl) 
     375            END DO 
     376         END DO 
     377 
     378         DO jj = 2, jpjm1                      !  Flux from i+1 to i IF u LT 0. 
     379            DO ji = fs_2, fs_jpim1 
     380               zbt  =       zbet(ji,jj) 
     381               zbt1 = 1.0 - zbet(ji,jj) 
     382               psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 
     383               zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 
     384               zalf1         = 1.0 - zalf 
     385               ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 
     386               ! 
     387               ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 
     388               psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) 
     389               psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 
     390                  &                                           + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) )    & 
     391                  &                                           + ( zalf1 - zalf ) * ztemp ) ) 
     392               psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  & 
     393                  &                                           + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) ) 
     394               psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) ) 
     395               psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) ) 
     396            END DO 
     397         END DO 
     398 
    314399      END DO 
    315400 
    316       !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                      
    317       psm (:,:)  = MAX( pcrh * e1e2t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 ) 
    318  
    319       !  Calculate fluxes and moments between boxes i<-->i+1               
    320       DO jj = 1, jpj                      !  Flux from i to i+1 WHEN u GT 0  
    321          DO ji = 1, jpi 
    322             zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 
    323             zalf         =  MAX( 0._wp, put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji,jj) 
    324             zalfq        =  zalf * zalf 
    325             zalf1        =  1.0 - zalf 
    326             zalf1q       =  zalf1 * zalf1 
    327             ! 
    328             zfm (ji,jj)  =  zalf  *   psm (ji,jj) 
    329             zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj) + zalf1 * ( psx(ji,jj) + (zalf1 - zalf) * psxx(ji,jj) )  ) 
    330             zfx (ji,jj)  =  zalfq * ( psx (ji,jj) + 3.0 * zalf1 * psxx(ji,jj) ) 
    331             zfxx(ji,jj)  =  zalf  *   psxx(ji,jj) * zalfq 
    332             zfy (ji,jj)  =  zalf  * ( psy (ji,jj) + zalf1 * psxy(ji,jj) ) 
    333             zfxy(ji,jj)  =  zalfq *   psxy(ji,jj) 
    334             zfyy(ji,jj)  =  zalf  *   psyy(ji,jj) 
    335  
    336             !  Readjust moments remaining in the box. 
    337             psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj) 
    338             ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj) 
    339             psx (ji,jj)  =  zalf1q * ( psx(ji,jj) - 3.0 * zalf * psxx(ji,jj) ) 
    340             psxx(ji,jj)  =  zalf1  * zalf1q * psxx(ji,jj) 
    341             psy (ji,jj)  =  psy (ji,jj) - zfy(ji,jj) 
    342             psyy(ji,jj)  =  psyy(ji,jj) - zfyy(ji,jj) 
    343             psxy(ji,jj)  =  zalf1q * psxy(ji,jj) 
    344          END DO 
    345       END DO 
    346  
    347       DO jj = 1, jpjm1                      !  Flux from i+1 to i when u LT 0. 
    348          DO ji = 1, fs_jpim1 
    349             zalf          = MAX( 0._wp, -put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji+1,jj)  
    350             zalg  (ji,jj) = zalf 
    351             zalfq         = zalf * zalf 
    352             zalf1         = 1.0 - zalf 
    353             zalg1 (ji,jj) = zalf1 
    354             zalf1q        = zalf1 * zalf1 
    355             zalg1q(ji,jj) = zalf1q 
    356             ! 
    357             zfm   (ji,jj) = zfm (ji,jj) + zalf  *   psm (ji+1,jj) 
    358             zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0 (ji+1,jj) - zalf1 * ( psx(ji+1,jj) - (zalf1 - zalf ) * psxx(ji+1,jj) ) ) 
    359             zfx   (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj) - 3.0 * zalf1 * psxx(ji+1,jj) ) 
    360             zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *   psxx(ji+1,jj) * zalfq 
    361             zfy   (ji,jj) = zfy (ji,jj) + zalf  * ( psy (ji+1,jj) - zalf1 * psxy(ji+1,jj) ) 
    362             zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *   psxy(ji+1,jj) 
    363             zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *   psyy(ji+1,jj) 
    364          END DO 
    365       END DO 
    366  
    367       DO jj = 2, jpjm1                     !  Readjust moments remaining in the box.  
    368          DO ji = fs_2, fs_jpim1 
    369             zbt  =       zbet(ji-1,jj) 
    370             zbt1 = 1.0 - zbet(ji-1,jj) 
    371             ! 
    372             psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji-1,jj) ) 
    373             ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji-1,jj) ) 
    374             psx (ji,jj) = zalg1q(ji-1,jj) * ( psx(ji,jj) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj) ) 
    375             psxx(ji,jj) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj) 
    376             psy (ji,jj) = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) - zfy (ji-1,jj) ) 
    377             psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) - zfyy(ji-1,jj) ) 
    378             psxy(ji,jj) = zalg1q(ji-1,jj) * psxy(ji,jj) 
    379          END DO 
    380       END DO 
    381  
    382       !   Put the temporary moments into appropriate neighboring boxes.     
    383       DO jj = 2, jpjm1                     !   Flux from i to i+1 IF u GT 0. 
    384          DO ji = fs_2, fs_jpim1 
    385             zbt  =       zbet(ji-1,jj) 
    386             zbt1 = 1.0 - zbet(ji-1,jj) 
    387             psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj) 
    388             zalf        = zbt * zfm(ji-1,jj) / psm(ji,jj) 
    389             zalf1       = 1.0 - zalf 
    390             ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji-1,jj) 
    391             ! 
    392             ps0 (ji,jj) = zbt * ( ps0(ji,jj) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj) 
    393             psx (ji,jj) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) + zbt1 * psx(ji,jj) 
    394             psxx(ji,jj) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj)                               & 
    395                &                + 5.0 * ( zalf * zalf1 * ( psx (ji,jj) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  )   & 
    396                &                                                + zbt1 * psxx(ji,jj) 
    397             psxy(ji,jj) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj)             & 
    398                &                + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj) ) )   & 
    399                &                                                + zbt1 * psxy(ji,jj) 
    400             psy (ji,jj) = zbt * ( psy (ji,jj) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj) 
    401             psyy(ji,jj) = zbt * ( psyy(ji,jj) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj) 
    402          END DO 
    403       END DO 
    404  
    405       DO jj = 2, jpjm1                     !  Flux from i+1 to i IF u LT 0. 
    406          DO ji = fs_2, fs_jpim1 
    407             zbt  =       zbet(ji,jj) 
    408             zbt1 = 1.0 - zbet(ji,jj) 
    409             psm(ji,jj)  = zbt * psm(ji,jj)  + zbt1 * ( psm(ji,jj) + zfm(ji,jj) ) 
    410             zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj) 
    411             zalf1       = 1.0 - zalf 
    412             ztemp       = - zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj) 
    413             ! 
    414             ps0(ji,jj)  = zbt * ps0 (ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) ) 
    415             psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) 
    416             psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( zalf * zalf * zfxx(ji,jj)  + zalf1 * zalf1 * psxx(ji,jj)  & 
    417                &                                      + 5.0 *( zalf * zalf1 * ( - psx(ji,jj) + zfx(ji,jj) )      & 
    418                &                                      + ( zalf1 - zalf ) * ztemp ) ) 
    419             psxy(ji,jj) = zbt * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)  & 
    420                &                                      + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj) )  ) 
    421             psy(ji,jj)  = zbt * psy (ji,jj)  + zbt1 * ( psy (ji,jj) + zfy (ji,jj) ) 
    422             psyy(ji,jj) = zbt * psyy(ji,jj)  + zbt1 * ( psyy(ji,jj) + zfyy(ji,jj) ) 
    423          END DO 
    424       END DO 
    425  
    426401      !-- Lateral boundary conditions 
    427       CALL lbc_lnk_multi( 'icedyn_adv_pra', psm , 'T',  1., ps0 , 'T',  1.   & 
    428          &              , psx , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes 
    429          &              , psxx, 'T',  1., psyy, 'T',  1.   & 
    430          &              , psxy, 'T',  1. ) 
    431  
    432       IF(ln_ctl) THEN 
    433          CALL prt_ctl(tab2d_1=psm  , clinfo1=' adv_x: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ') 
    434          CALL prt_ctl(tab2d_1=psx  , clinfo1=' adv_x: psx  :', tab2d_2=psxx, clinfo2=' psxx : ') 
    435          CALL prt_ctl(tab2d_1=psy  , clinfo1=' adv_x: psy  :', tab2d_2=psyy, clinfo2=' psyy : ') 
    436          CALL prt_ctl(tab2d_1=psxy , clinfo1=' adv_x: psxy :') 
    437       ENDIF 
     402      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1., ps0 , 'T',  1.   & 
     403         &                                , psx             , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes 
     404         &                                , psxx            , 'T',  1., psyy, 'T',  1. , psxy, 'T',  1. ) 
    438405      ! 
    439406   END SUBROUTINE adv_x 
    440407 
    441408 
    442    SUBROUTINE adv_y( pdf, pvt , pcrh, psm , ps0 ,   & 
     409   SUBROUTINE adv_y( pdt, pvt , pcrh, psm , ps0 ,   & 
    443410      &              psx, psxx, psy , psyy, psxy ) 
    444411      !!--------------------------------------------------------------------- 
     
    448415      !!                variable on y axis 
    449416      !!--------------------------------------------------------------------- 
    450       REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step 
    451       REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0) 
    452       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s] 
    453       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area 
    454       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected 
    455       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments  
    456       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments 
     417      REAL(wp)                  , INTENT(in   ) ::   pdt                ! time step 
     418      REAL(wp)                  , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0) 
     419      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s] 
     420      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psm                ! area 
     421      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ps0                ! field to be advected 
     422      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psx , psy          ! 1st moments  
     423      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments 
    457424      !! 
    458       INTEGER  ::   ji, jj                               ! dummy loop indices 
    459       REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp         ! temporary scalars 
     425      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices 
     426      REAL(wp) ::   zs1max, zslpmax, ztemp               ! temporary scalars 
    460427      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         - 
    461428      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         - 
     
    464431      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      - 
    465432      !--------------------------------------------------------------------- 
    466  
    467       ! Limitation of moments. 
    468  
    469       zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection. 
    470  
    471       DO jj = 1, jpj 
    472          DO ji = 1, jpi 
    473             zslpmax = MAX( 0._wp, ps0(ji,jj) ) 
    474             zs1max  = 1.5 * zslpmax 
    475             zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) ) 
    476             zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   & 
    477                &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  ) 
    478             rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
    479             ! 
    480             ps0 (ji,jj) = zslpmax   
    481             psx (ji,jj) = psx (ji,jj) * rswitch 
    482             psxx(ji,jj) = psxx(ji,jj) * rswitch 
    483             psy (ji,jj) = zs1new * rswitch 
    484             psyy(ji,jj) = zs2new * rswitch 
    485             psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch 
    486          END DO 
     433      ! 
     434      jcat = SIZE( ps0 , 3 )   ! size of input arrays 
     435      !       
     436      DO jl = 1, jcat   ! loop on categories 
     437         ! 
     438         ! Limitation of moments. 
     439         DO jj = 1, jpj 
     440            DO ji = fs_2, fs_jpim1 
     441               !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 
     442               psm(ji,jj,jl) = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20  ) 
     443               ! 
     444               zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 
     445               zs1max  = 1.5 * zslpmax 
     446               zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) ) 
     447               zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   & 
     448                  &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) )  ) 
     449               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
     450               ! 
     451               ps0 (ji,jj,jl) = zslpmax   
     452               psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch 
     453               psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch 
     454               psy (ji,jj,jl) = zs1new         * rswitch 
     455               psyy(ji,jj,jl) = zs2new         * rswitch 
     456               psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 
     457            END DO 
     458         END DO 
     459  
     460         !  Calculate fluxes and moments between boxes j<-->j+1               
     461         DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0    
     462            DO ji = fs_2, fs_jpim1 
     463               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 
     464               zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt * e1v(ji,jj) / psm(ji,jj,jl) 
     465               zalfq        =  zalf * zalf 
     466               zalf1        =  1.0 - zalf 
     467               zalf1q       =  zalf1 * zalf1 
     468               ! 
     469               zfm (ji,jj)  =  zalf  * psm(ji,jj,jl) 
     470               zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl)  + (zalf1-zalf) * psyy(ji,jj,jl) ) )  
     471               zfy (ji,jj)  =  zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) ) 
     472               zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj,jl) 
     473               zfx (ji,jj)  =  zalf  * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 
     474               zfxy(ji,jj)  =  zalfq * psxy(ji,jj,jl) 
     475               zfxx(ji,jj)  =  zalf  * psxx(ji,jj,jl) 
     476               ! 
     477               !  Readjust moments remaining in the box. 
     478               psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj) 
     479               ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj) 
     480               psy (ji,jj,jl)  =  zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) ) 
     481               psyy(ji,jj,jl)  =  zalf1 * zalf1q * psyy(ji,jj,jl) 
     482               psx (ji,jj,jl)  =  psx (ji,jj,jl) - zfx(ji,jj) 
     483               psxx(ji,jj,jl)  =  psxx(ji,jj,jl) - zfxx(ji,jj) 
     484               psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl) 
     485            END DO 
     486         END DO 
     487         ! 
     488         DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0. 
     489            DO ji = fs_2, fs_jpim1 
     490               zalf          = ( MAX(0._wp, -pvt(ji,jj) ) * pdt * e1v(ji,jj) ) / psm(ji,jj+1,jl)  
     491               zalg  (ji,jj) = zalf 
     492               zalfq         = zalf * zalf 
     493               zalf1         = 1.0 - zalf 
     494               zalg1 (ji,jj) = zalf1 
     495               zalf1q        = zalf1 * zalf1 
     496               zalg1q(ji,jj) = zalf1q 
     497               ! 
     498               zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji,jj+1,jl) 
     499               zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji,jj+1,jl) & 
     500                  &                                   - zalf1 * (psy(ji,jj+1,jl) - (zalf1 - zalf ) * psyy(ji,jj+1,jl) ) ) 
     501               zfy   (ji,jj) = zfy (ji,jj) + zalfq * (  psy (ji,jj+1,jl) - 3.0 * zalf1 * psyy(ji,jj+1,jl) ) 
     502               zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji,jj+1,jl) * zalfq 
     503               zfx   (ji,jj) = zfx (ji,jj) + zalf  * (  psx (ji,jj+1,jl) - zalf1 * psxy(ji,jj+1,jl) ) 
     504               zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji,jj+1,jl) 
     505               zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji,jj+1,jl) 
     506            END DO 
     507         END DO 
     508 
     509         !  Readjust moments remaining in the box.  
     510         DO jj = 2, jpjm1 
     511            DO ji = fs_2, fs_jpim1 
     512               zbt  =         zbet(ji,jj-1) 
     513               zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
     514               ! 
     515               psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) ) 
     516               ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) ) 
     517               psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) ) 
     518               psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl) 
     519               psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) ) 
     520               psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) ) 
     521               psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl) 
     522            END DO 
     523         END DO 
     524 
     525         !   Put the temporary moments into appropriate neighboring boxes.     
     526         DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0. 
     527            DO ji = fs_2, fs_jpim1 
     528               zbt  =       zbet(ji,jj-1) 
     529               zbt1 = 1.0 - zbet(ji,jj-1) 
     530               psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl)  
     531               zalf          = zbt * zfm(ji,jj-1) / psm(ji,jj,jl)  
     532               zalf1         = 1.0 - zalf 
     533               ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1) 
     534               ! 
     535               ps0(ji,jj,jl)  =   zbt  * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl) 
     536               psy(ji,jj,jl)  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )  & 
     537                  &             + zbt1 * psy(ji,jj,jl)   
     538               psyy(ji,jj,jl) =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl)                           & 
     539                  &                      + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &  
     540                  &             + zbt1 * psyy(ji,jj,jl) 
     541               psxy(ji,jj,jl) =   zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl)            & 
     542                  &                      + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) )  & 
     543                  &             + zbt1 * psxy(ji,jj,jl) 
     544               psx (ji,jj,jl) =   zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl) 
     545               psxx(ji,jj,jl) =   zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl) 
     546            END DO 
     547         END DO 
     548 
     549         DO jj = 2, jpjm1                      !  Flux from j+1 to j IF v LT 0. 
     550            DO ji = fs_2, fs_jpim1 
     551               zbt  =       zbet(ji,jj) 
     552               zbt1 = 1.0 - zbet(ji,jj) 
     553               psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 
     554               zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 
     555               zalf1         = 1.0 - zalf 
     556               ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 
     557               ! 
     558               ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * (  ps0(ji,jj,jl) + zf0(ji,jj) ) 
     559               psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) 
     560               psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) & 
     561                  &                                            + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) )    & 
     562                  &                                            + ( zalf1 - zalf ) * ztemp ) ) 
     563               psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  & 
     564                  &                                            + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) ) 
     565               psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) ) 
     566               psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) ) 
     567            END DO 
     568         END DO 
     569 
    487570      END DO 
    488571 
    489       !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 
    490       psm(:,:)  = MAX(  pcrh * e1e2t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20  ) 
    491  
    492       !  Calculate fluxes and moments between boxes j<-->j+1               
    493       DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0    
    494          DO ji = 1, jpi 
    495             zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 
    496             zalf         =  MAX( 0._wp, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj) 
    497             zalfq        =  zalf * zalf 
    498             zalf1        =  1.0 - zalf 
    499             zalf1q       =  zalf1 * zalf1 
    500             ! 
    501             zfm (ji,jj)  =  zalf  * psm(ji,jj) 
    502             zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj)  + (zalf1-zalf) * psyy(ji,jj)  ) )  
    503             zfy (ji,jj)  =  zalfq *( psy(ji,jj) + 3.0*zalf1*psyy(ji,jj) ) 
    504             zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj) 
    505             zfx (ji,jj)  =  zalf  * ( psx(ji,jj) + zalf1 * psxy(ji,jj) ) 
    506             zfxy(ji,jj)  =  zalfq * psxy(ji,jj) 
    507             zfxx(ji,jj)  =  zalf  * psxx(ji,jj) 
    508             ! 
    509             !  Readjust moments remaining in the box. 
    510             psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj) 
    511             ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj) 
    512             psy (ji,jj)  =  zalf1q * ( psy(ji,jj) -3.0 * zalf * psyy(ji,jj) ) 
    513             psyy(ji,jj)  =  zalf1 * zalf1q * psyy(ji,jj) 
    514             psx (ji,jj)  =  psx (ji,jj) - zfx(ji,jj) 
    515             psxx(ji,jj)  =  psxx(ji,jj) - zfxx(ji,jj) 
    516             psxy(ji,jj)  =  zalf1q * psxy(ji,jj) 
    517          END DO 
    518       END DO 
    519       ! 
    520       DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0. 
    521          DO ji = 1, jpi 
    522             zalf          = ( MAX(0._wp, -pvt(ji,jj) ) * zrdt * e1v(ji,jj) ) / psm(ji,jj+1)  
    523             zalg  (ji,jj) = zalf 
    524             zalfq         = zalf * zalf 
    525             zalf1         = 1.0 - zalf 
    526             zalg1 (ji,jj) = zalf1 
    527             zalf1q        = zalf1 * zalf1 
    528             zalg1q(ji,jj) = zalf1q 
    529             ! 
    530             zfm   (ji,jj) = zfm (ji,jj) + zalf  *   psm (ji,jj+1) 
    531             zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0 (ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) ) 
    532             zfy   (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) ) 
    533             zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *   psyy(ji,jj+1) * zalfq 
    534             zfx   (ji,jj) = zfx (ji,jj) + zalf  * ( psx (ji,jj+1) - zalf1 * psxy(ji,jj+1) ) 
    535             zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *   psxy(ji,jj+1) 
    536             zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *   psxx(ji,jj+1) 
    537          END DO 
    538       END DO 
    539  
    540       !  Readjust moments remaining in the box.  
    541       DO jj = 2, jpj 
    542          DO ji = 1, jpi 
    543             zbt  =         zbet(ji,jj-1) 
    544             zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
    545             ! 
    546             psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) ) 
    547             ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) ) 
    548             psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) ) 
    549             psyy(ji,jj) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj) 
    550             psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) ) 
    551             psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) ) 
    552             psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj) 
    553          END DO 
    554       END DO 
    555  
    556       !   Put the temporary moments into appropriate neighboring boxes.     
    557       DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0. 
    558          DO ji = 1, jpi 
    559             zbt  =         zbet(ji,jj-1) 
    560             zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
    561             psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj)  
    562             zalf        = zbt * zfm(ji,jj-1) / psm(ji,jj)  
    563             zalf1       = 1.0 - zalf 
    564             ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1) 
    565             ! 
    566             ps0(ji,jj)  = zbt  * ( ps0(ji,jj) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj) 
    567             psy(ji,jj)  = zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp )   & 
    568                &                                               + zbt1 * psy(ji,jj)   
    569             psyy(ji,jj) = zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj)                             & 
    570                &                 + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )   &  
    571                &                                               + zbt1 * psyy(ji,jj) 
    572             psxy(ji,jj) = zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj)               & 
    573                &                  + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) )  )   & 
    574                &                                                + zbt1 * psxy(ji,jj) 
    575             psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj) 
    576             psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj) 
    577          END DO 
    578       END DO 
    579  
    580       DO jj = 2, jpjm1                   !  Flux from j+1 to j IF v LT 0. 
    581          DO ji = 1, jpi 
    582             zbt  =         zbet(ji,jj) 
    583             zbt1 = ( 1.0 - zbet(ji,jj) ) 
    584             psm(ji,jj)  = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) ) 
    585             zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj) 
    586             zalf1       = 1.0 - zalf 
    587             ztemp       = - zalf * ps0 (ji,jj) + zalf1 * zf0(ji,jj) 
    588             ps0 (ji,jj) =   zbt  * ps0 (ji,jj) + zbt1  * ( ps0(ji,jj) + zf0(ji,jj) ) 
    589             psy (ji,jj) =   zbt  * psy (ji,jj) + zbt1  * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp ) 
    590             psyy(ji,jj) =   zbt  * psyy(ji,jj) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj)   & 
    591                &                                         + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) )          & 
    592                &                                         + ( zalf1 - zalf ) * ztemp )                                ) 
    593             psxy(ji,jj) =   zbt  * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)       & 
    594                &                                         + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) )  ) 
    595             psx (ji,jj) =   zbt  * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) ) 
    596             psxx(ji,jj) =   zbt  * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) ) 
    597          END DO 
    598       END DO 
    599  
    600572      !-- Lateral boundary conditions 
    601       CALL lbc_lnk_multi( 'icedyn_adv_pra', psm , 'T',  1.,  ps0 , 'T',  1.   & 
    602          &              , psx , 'T', -1.,  psy , 'T', -1.   &   ! caution gradient ==> the sign changes 
    603          &              , psxx, 'T',  1.,  psyy, 'T',  1.   & 
    604          &              , psxy, 'T',  1. ) 
    605  
    606       IF(ln_ctl) THEN 
    607          CALL prt_ctl(tab2d_1=psm  , clinfo1=' adv_y: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ') 
    608          CALL prt_ctl(tab2d_1=psx  , clinfo1=' adv_y: psx  :', tab2d_2=psxx, clinfo2=' psxx : ') 
    609          CALL prt_ctl(tab2d_1=psy  , clinfo1=' adv_y: psy  :', tab2d_2=psyy, clinfo2=' psyy : ') 
    610          CALL prt_ctl(tab2d_1=psxy , clinfo1=' adv_y: psxy :') 
    611       ENDIF 
     573      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1., ps0 , 'T',  1.   & 
     574         &                                , psx             , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes 
     575         &                                , psxx            , 'T',  1., psyy, 'T',  1. , psxy, 'T',  1. ) 
    612576      ! 
    613577   END SUBROUTINE adv_y 
     
    624588      ! 
    625589      !                             !* allocate prather fields 
    626       ALLOCATE( sxopw(jpi,jpj)     , syopw(jpi,jpj)     , sxxopw(jpi,jpj)     , syyopw(jpi,jpj)     , sxyopw(jpi,jpj)     ,   & 
     590      ALLOCATE( sxopw(jpi,jpj,1)   , syopw(jpi,jpj,1)   , sxxopw(jpi,jpj,1)   , syyopw(jpi,jpj,1)   , sxyopw(jpi,jpj,1)   ,   & 
    627591         &      sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) ,   & 
    628592         &      sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) ,   & 
     
    652616      !!                   ***  ROUTINE adv_pra_rst  *** 
    653617      !!                      
    654       !! ** Purpose :   Read or write RHG file in restart file 
     618      !! ** Purpose :   Read or write file in restart file 
    655619      !! 
    656620      !! ** Method  :   use of IOM library 
  • NEMO/trunk/src/ICE/icedyn_adv_umx.F90

    r10945 r11612  
    319319         ! 
    320320         !== Ice age ==! 
    321          IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN 
    322             zamsk = 1._wp 
    323             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, & 
    324                &                                      poa_i, poa_i ) 
    325          ENDIF 
     321         zamsk = 1._wp 
     322         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, & 
     323            &                                      poa_i, poa_i ) 
    326324         ! 
    327325         !== melt ponds ==! 
Note: See TracChangeset for help on using the changeset viewer.